TRCA-OSAP Stream Habitat Analysis

Introduction

Background

The ability of streams and rivers to support a healthy and diverse biota depends on the availability of suitable habitat. This includes channel structure (e.g., cover, substrate) and channel processes (e.g., hydrology and sediment transport). Understanding these habitat components of a river system is crucial as it influences habitat availability, water quality, and overall ecosystem health. In the context of Ontario Stream Assessment Protocol (OSAP) and the Toronto and Region Conservation Authority (TRCA), in-stream habitat monitoring can provide valuable insights into the health of aquatic ecosystems, particularly for fish and benthic invertebrates, which are often sensitive to changes in their physical environment. By examining the relationships between habitat conditions and biological responses, we can better understand the impacts of natural and anthropogenic disturbances on aquatic ecosystems to guide land-use planning and restoration efforts.

Project Overview

This study aims to understanding the in-stream habitat data collected through TRCA’s long-term monitoring programs. It includes compilation and summarization of the data, as well as exploring spatial and temporal patterns. The results will allow TRCA and OSAP to assess the variability in habitat across the region and how they relate to fish and benthic invertebrate communities to better inform monitoring and management decisions.

TRCA Monitoring Locations

Methods

Data and models

Habitat and fish community data were downloaded from the Flowing Waters Information System (FWIS) Channel Structure Summary table. Watershed outlines were downloaded from TRCA open data portal. Supplemental data was joined from V3 of the Ontario Aquatic Ecosystem Classification (AEC). The Upstream catchment stress index was used to generate 3 Stress Classes corresponding to low (20-80), medium (80-140), and high (140-200) anthropogenic impacts. The Upstream Catchment Area (UCA) was used to generate 3 UCA Classes corresponding to small (0-20), intermediate (20-80), and large (>80) stream sizes.

Generalized Linear Mixed Effects Models (GLMMs; using the GLMMtmb R package) were used to describe spatial and temporal patterns in in-stream habitat endpoints over time. Distributional families were selected to match the distributions and data generating processes that underlie the endpoints, for example substrate sizes tend to vary on a log-normal distribution, whereas percentages and cover endpoints generally follow a beta distribution.

Some endpoints contained multiple categorical values within continuously measured values (i.e., sand, silt, clay categories for fines within the Substrate Sizes endpoints). These values were all replaced by 0’s and modeled with the zero-inflation portion of the model. See the R Code for full list of data transformations. Endpoints with more than 75% of observations at 0 tended to have worse preliminary model fits and were excluded from the final analysis.

Each endpoint was modeled with the following formula: \[ Endpoint = ln(Reach Slope+1) + Year_n*UCA*Stress + (Year_n|SiteCode) + (1|Year_n) \]

where \(Reach Slope\) is the percent slope of the reach from the AEC; \(UCA\), and \(Stress\) are fixed effect factor terms; \(Year_n\) represent temporal slopes and the \((Year_n|SiteCode)\) is a nested random-effects term which fits a temporal slope to each monitoring location nested within a UCA class, and Stress Class factor. The \(Year_n*UCA*Stress\) term is fit as three-way interaction terms fitting temporal slopes to the two factor variables together and individually. The \((1|Year_n)\) term fits a random intercept to each monitoring year. This allows for global effects that can influence all sites in a particular year such as a potential crew bias or severe weather event.

Some endpoints are also fit with a zero-inflation model (probability of a zero: \(E_0\)). Zero-inflation was a simplified model without the three-way interaction \(Year_n*UCA*Stress\) term:

\[ E_0 = ln(Reach Slope+1) + Year_n + UCA + Stress + (Year_n|SiteCode) + (1|Year_n) \]

The following table summarizes the distributional family used for each endpoint.

	Distribution
Cover
Wood Cover - Unembedded	beta
Wood Cover - Total	beta
Round Rock Cover - Unembedded	beta
Round Rock Cover - Total	beta
Flat Rock Cover - Unembedded	beta
Flat Rock Cover - Total	beta
Percent
Filamentous Algae Percent	beta
Non-Filamentous Algae Percent	beta
Moss Percent	beta
Indices
Sediment Sorting Index	gamma
Bank Stability - Mean	beta
Habitat Homogeneity	beta
Substrate Sizes
D50	lognormal
D85	lognormal
D50Max	lognormal
Morphology
Width - Mean	lognormal
Depth - Mean	lognormal
Width:Depth Ratio	gamma
Usable Area	gamma

After an initial round of model fitting, models with convergence errors were subject to an outlier removal method. 10,000 simulated residuals were generated from the fitted models, and observed residuals were compared to simulated residuals using the DHARMa R package. Observed residuals that were outside the 2.5th and 97.5th percentiles of simulated residuals were considered outliers. These values were removed from the model training data, and models were refit to the reduced dataset. Outlier values were excluded from model fitting, but included in final evaluations.

R Code

This section contains all the R code to reproduce the current results and analysis.

0 FWIS Source

if( !require( "httr" ) ) {
  install.packages( "httr" )
}
library( httr )
if( !require( "jsonlite" ) ) {
  install.packages( "jsonlite" )
}
library( jsonlite )
if( !require( "base64enc" ) ) {
  install.packages( "base64enc" )
}
library( base64enc )

# widePost.r -- a collection of functions that enable access to WIDE server data using R
#
# Usage: source( "widePost.r" )  ...from R console/script, with current directory set to folder with widePost.r
#
# to enable some debugging, such as the HTTP request/response headers and JSON request, enter the following:
#  options( verbose = 1 )
# to disable subsequent debugging, use:
#  options( verbose = 0 )

.WIDE <- new.env()
.WIDE$clientVersionDate <- "2016-06-27-1400"
.WIDE$clientVersion <- c( "1", "0" )
# list of servers that the user has logged in to...
.WIDE$connectedServers <- list()

# wideLogin -- connect to a WIDE data API server
#  serverURL: the URL of the server data API script
#  username/password: login credentials
# returns false in the event of an error, or a serverHandle (integer > 0) value to be passed to subsequent widePost functions

wideLogin<-function( serverURL, username, password )
{
  # application.php requires parameters "Username", "Password" and "cls_logon_submit" as follows...
  postBody<-list( Username = username, Password = password, cls_logon_submit = 1 )
  loginResponse <- httr::POST( serverURL, body = postBody, encode = "form" )
  loginContent <- content( loginResponse, "text" )
  
  # error handling
  if( !jsonlite::validate( loginContent ) ) {
    # php syntax errors do not get encoded as JSON, so check for beginning and ending curly braces...
    if( ( substr( loginContent, 1, 1 ) != "{" ) &&
        ( substr( loginContent, nchar( loginContent ), nchar( loginContent ) ) != "}" ) ) {
      print( paste( " php syntax error: ", loginContent ) )
    } else {
      print( paste( " invalid JSON: ", loginContent ) )
    }
    return( FALSE )
  }
  
  jsonData <- fromJSON( loginContent )
  
  if( "error" %in% names( jsonData ) ) {
    print( jsonData )
    return( FALSE )
  }
  
  # successful login
  print( jsonData$status )
  
  # check for version differences
  if( .WIDE$clientVersion[ 1 ] == jsonData$version[ 1 ] ) {
    if( .WIDE$clientVersion[ 2 ] == jsonData$version[ 2 ] ) {
      print( paste( "(0) WIDE interface ready - client/server version", paste( .WIDE$clientVersion, collapse = "." ),
                    " (", .WIDE$clientVersionDate, ")" ) )
      versionStatus <- 0
    } else if( .WIDE$clientVersion[ 2 ] > jsonData$version[ 2 ] ) {
      versionStatus <- 1
      print( paste0( "(1) WIDE interface ready - server version ", jsonData$version[ 1 ], ".", jsonData$version[ 2 ],
                     " (", jsonData$versionDate, "), client: ", paste( .WIDE$clientVersion, collapse = "." ), " (",
                     .WIDE$clientVersionDate, ")" ) )
    } else {
      versionStatus <- -1
      cat( paste0( "(-1) WIDE interface ready - server version ", jsonData$version[ 1 ], ".", jsonData$version[ 2 ],
                   " (", jsonData$versionDate, "), client: ", paste( .WIDE$clientVersion, collapse = "." ), " (",
                   .WIDE$clientVersionDate, ")\nPlease consider downloading a new client\n" ) )
    }
  } else if( .WIDE$clientVersion[ 1 ] > jsonData$version[ 1 ] ) {
    versionStatus <- 1
    print( paste0( "(1) WIDE interface ready - server version ", jsonData$version[ 1 ], ".", jsonData$version[ 2 ],
                   " (", jsonData$versionDate, "), client: ", paste( .WIDE$clientVersion, collapse = "." ), " (", .WIDE$clientVersionDate, ")" ) )
  } else {
    versionStatus <- -1
    print( paste0( "(-1) FATAL ERROR: client is out of date - server version ", jsonData$version[ 1 ], ".", jsonData$version[ 2 ],
                   " (", jsonData$versionDate, "), client: ", paste( .WIDE$clientVersion, collapse = "." ), " (", .WIDE$clientVersionDate,
                   ")\nPlease consider downloading a new client\n" ) )
  }
  
  # create and store server handle to be returned to the user
  if( versionStatus == -1 ) {
    return( -1 )
  }
  .WIDE$connectedServers <- append( .WIDE$connectedServers, serverURL )
  return( length( .WIDE$connectedServers ) )
}

# wideDataSelect -- retrieve one or more records from a database table
#  e.g.,
#   wideDataSelect(serverHandle,"lkpStreamCode",list("StreamName"))
#   wideDataSelect(serverHandle,"lkpStreamCode",list("StreamName","StreamCode"),orderBy=list("StreamName"))
#   wideDataSelect(serverHandle,"lkpStreamCode",list("StreamName","StreamCode"),orderBy=list("StreamName ASC"))
#   wideDataSelect(serverHandle,"lkpStreamCode",list("StreamName","StreamCode"),orderBy=list("StreamName DESC"))
#   wideDataSelect(serverHandle,"lkpStreamCode",list("StreamName","StreamCode"),orderBy=list('"1" DESC'))
#   wideDataSelect(serverHandle,"lkpStreamCode",list("StreamName","StreamCode"),orderBy=list('1 DESC'))
#   wideDataSelect(serverHandle,"lkpStreamCode",list("StreamName","StreamCode"),orderBy=list(1))
#   wideDataSelect(serverHandle,"lkpStreamCode",list("StreamName","StreamCode"),list(StreamName=c("like","%Black%")))
#   wideDataSelect(serverHandle,"lkpStreamCode",list("StreamName","StreamCode"),list(StreamName=c("like","%Black%"),StreamCode=c("like","bl%")),list(1),"OR")
# returns a data.frame if successful; otherwise, returns NULL

wideDataSelect<-function( serverHandle, table, colsSelected, where = NULL, orderBy = NULL, logical = "AND", geoConstraint = 0 )
{
  if( getOption( "wideFunCall", 0 ) ) {
    print( call( "wideDataSelect", serverHandle = serverHandle, table = table, colsSelected = colsSelected,
                 where = where, orderBy = orderBy, logical = logical, geoConstraint = geoConstraint ) )
  }
  if( ( serverHandle < 1 ) || ( serverHandle > length( .WIDE$connectedServers ) ) ) {
    stop( paste0( "Invalid server handle, must be one of: [", seq_len( length( .WIDE$connectedServers ) ), "]" ) )
  }
  
  # Note that the "xxx=" prefix in each "element" of the jsonPost list must be specified (i.e., it can't be left out
  #  of the jsonPost list() below)...
  # e.g., jsonPost<-list(table="lkpStreamCode",action="select",colsSelected=list("StreamName"))
  jsonPost<-list( table = table, action = "select", colsSelected = colsSelected, colsWhere = where, orderbyCols = orderBy, logical = logical, geoConstraint = geoConstraint )
  
  if( getOption( "verbose" ) ) {
    print( jsonPost )
    postResponse<-httr::POST( .WIDE$connectedServers[[serverHandle]], body = jsonPost, encode = "json", verbose() )
  } else {
    postResponse<-httr::POST( .WIDE$connectedServers[[serverHandle]], body = jsonPost, encode = "json" )
  }
  # for some circumstances (e.g., application development, testing, etc.) it may be useful to save the basic postResponse content
  #  to a file...
  # write( paste( "Content of postResponse ", content( postResponse, "text" ) ), "c:/temp/routput.txt" )
  
  postContent<-content( postResponse, "text" )
  
  if( !jsonlite::validate( postContent ) ) {
    if( ( substr( postContent, 1, 1 ) != "{" ) &&
        ( substr( postContent, nchar( postContent ), nchar( postContent ) ) != "}" ) ) {
      print( paste( " php syntax error: ", postContent ) )
    } else {
      print( paste( " invalid json: ", postContent ) )
    }
    return( NULL )
  }
  
  jsonData <- fromJSON( postContent )
  
  if( "error" %in% names( jsonData ) ) {
    print( jsonData )
    return( NULL )
  }
  
  # find all of the binary columns which must be decoded
  binaryCols <- Filter( function( colName ) { jsonData[[ 1, colName ]] == 1 }, names( jsonData ) )
  
  # removing a row from a dataframe with one column ends up converting that data frame into a vector;
  # to prevent this from happening, use drop=FALSE
  jsonData <- jsonData[ -1, , drop = FALSE ]
  
  # if the data frame isn't empty after removing the row, renumber all of the rows
  if( nrow( jsonData ) ) {
    row.names( jsonData ) <- 1:nrow( jsonData )
  }
  
  # if there are any binary columns to decode, do that now... (base64decode cannot NA, so those must be explicitly
  #  returned back as NA)
  # NA_character_ is returned because that is the type of NA value used in the shinyapp DataTable
  if( length( binaryCols ) ) {
    jsonData[[ binaryCols ]] <- sapply( jsonData[[ binaryCols ]], function( i ) {
      if( is.na( i ) ) {
        return( NA_character_ )
      }
      # for some reason, without the paste (... collapse), the as.character(...) ends up as type list
      #  in the final jsonData dataframe
      return( paste0( "0x", as.character( base64decode( i ) ), collapse = "" ) )
    } )
  }
  
  return( jsonData )
}

# wideDataSpatialConstraint -- define a geospatial constraint for subsequent queries with
#   spatial constraints enabled
#  e.g.,

# return TRUE if successful; otherwise, return FALSE
wideDataSpatialConstraint <- function( serverHandle, table = NULL, column = NULL, where = NULL, comparisonResult = NULL, logical="AND") {
  
  if( getOption( "wideFunCall", 0 ) ) {
    print( call( "wideDataSpatialConstraint", serverHandle = serverHandle, table = table, column = column, where = where,
                 comparisonResult = comparisonResult, logical = logical ) )
  }
  if( ( serverHandle < 1 ) || ( serverHandle > length( .WIDE$connectedServers ) ) ) {
    stop( paste0( "Invalid server handle, must be one of: ", seq_len( length( .WIDE$connectedServers ) ) ) )
  }
  
  jsonPost<-list( table = table, action = "geospatialConstraint", column = column, colsWhere = where,
                  comparisonResult = comparisonResult, logical = logical )
  
  if( getOption( "verbose" ) ) {
    print( jsonPost )
    postResponse<-httr::POST( .WIDE$connectedServers[[serverHandle]], body = jsonPost, encode = "json", verbose() )
  } else {
    postResponse<-httr::POST( .WIDE$connectedServers[[serverHandle]], body = jsonPost, encode = "json" )
  }
  
  postContent<-content( postResponse, "text" )
  
  if( !jsonlite::validate( postContent ) ) {
    if( ( substr( postContent, 1, 1 ) != "{" ) &&
        ( substr( postContent, nchar( postContent ), nchar( postContent ) ) != "}" ) ) {
      print( paste( " php syntax error: ", postContent ) )
    } else {
      print( paste( " invalid json: ", postContent ) )
    }
    return( FALSE )
  }
  
  jsonData <- fromJSON( postContent )
  print( jsonData )
  
  if( "error" %in% names( jsonData ) ) {
    return( FALSE )
  }
  return( TRUE )
}

# wideDataUpdate -- update one or more records of a database table
#  e.g.,
# wideDataUpdate(serverHandle = 1L, table = "lkpStreamCode", cols = list(StreamName = "Teststream1"), where =list(StreamCode = c("=", "teststr")))
# return TRUE if successful; otherwise, return FALSE

wideDataUpdate<-function( serverHandle, table, cols, where, logical = "AND" )
{
  if( getOption( "wideFunCall", 0 ) ) {
    print( call( "wideDataUpdate", serverHandle = serverHandle, table = table, cols = cols, where = where, logical = logical ) )
  }
  if( ( serverHandle < 1 ) || ( serverHandle > length( .WIDE$connectedServers ) ) ) {
    stop( paste0( "Invalid server handle, must be one of: ", seq_len( length( .WIDE$connectedServers ) ) ) )
  }
  
  # Note that the "xxx=" prefix in each "element" of the jsonPost list must be specified (i.e., it can't be left out of the jsonPost list() below)...
  jsonPost<-list( table = table, action = "update", cols = cols, colsWhere = where, logical = logical )
  
  if( getOption( "verbose" ) ) {
    print( jsonPost )
    postResponse<-httr::POST( .WIDE$connectedServers[[serverHandle]], body = jsonPost, encode = "json", verbose() )
  } else {
    postResponse<-httr::POST( .WIDE$connectedServers[[serverHandle]], body = jsonPost, encode = "json" )
  }
  
  postContent<-content( postResponse, "text" )
  
  if( !jsonlite::validate( postContent ) ) {
    if( ( substr( postContent, 1, 1 ) != "{" ) &&
        ( substr( postContent, nchar( postContent ), nchar( postContent ) ) != "}" ) ) {
      print( paste( " php syntax error: ", postContent ) )
    } else {
      print( paste( " invalid json: ", postContent ) )
    }
    return( FALSE )
  }
  
  jsonData <- fromJSON( postContent )
  print( jsonData )
  
  if( "error" %in% names( jsonData ) ) {
    return( FALSE )
  }
  
  return( TRUE )
}

# return TRUE if successful; otherwise, return FALSE

wideDataInsert<-function( serverHandle, table, cols )
{
  if( getOption( "wideFunCall", 0 ) ) {
    print( call( "wideDataInsert", serverHandle = serverHandle, table = table, cols = cols ) )
  }
  if( ( serverHandle < 1 ) || ( serverHandle > length( .WIDE$connectedServers ) ) ) {
    stop( paste0( "Invalid server handle, must be one of: ", seq_len( length( .WIDE$connectedServers ) ) ) )
  }
  
  jsonPost<-list( table = table, action = "insert", cols = cols )
  
  if( getOption( "verbose" ) ) {
    print( jsonPost )
    postResponse<-httr::POST( .WIDE$connectedServers[[ serverHandle ]], body = jsonPost, encode = "json", verbose() )
  } else {
    postResponse<-httr::POST( .WIDE$connectedServers[[ serverHandle ]], body = jsonPost, encode = "json" )
  }
  
  postContent<-content( postResponse, "text" )
  
  if( !jsonlite::validate( postContent ) ) {
    if( ( substr( postContent, 1, 1 ) != "{" ) &&
        ( substr( postContent, nchar( postContent ), nchar( postContent ) ) != "}" ) ) {
      print( paste( " php syntax error: ", postContent ) )
    } else {
      print( paste( " invalid JSON: ", postContent ) )
    }
    return( FALSE )
  }
  
  jsonData <- fromJSON( postContent )
  print( jsonData )
  
  if( "error" %in% names( jsonData ) ) {
    return( FALSE )
  }
  return( TRUE )
}

# return TRUE if successful; otherwise, return FALSE

wideDataDelete<-function( serverHandle, table, where, logical = "AND" )
{
  if( getOption( "wideFunCall", 0 ) ) {
    print( call( "wideDataDelete", serverHandle = serverHandle, table = table, where = where, logical = logical ) )
  }
  if( ( serverHandle < 1 ) || ( serverHandle > length( .WIDE$connectedServers ) ) ) {
    stop( paste0( "Invalid server handle, must be one of: ", seq_len( length( .WIDE$connectedServers ) ) ) )
  }
  
  jsonPost<-list( table = table, action = "delete", colsWhere = where, logical = logical )
  
  if( getOption( "verbose" ) ) {
    print( jsonPost )
    postResponse<-httr::POST( .WIDE$connectedServers[[ serverHandle ]], body = jsonPost, encode = "json", verbose() )
  } else {
    postResponse<-httr::POST( .WIDE$connectedServers[[ serverHandle ]], body = jsonPost, encode = "json" )
  }
  
  postContent<-content( postResponse, "text" )
  if( !jsonlite::validate( postContent ) ) {
    if( ( substr( postContent, 1, 1 ) != "{" ) &&
        ( substr( postContent, nchar( postContent ), nchar( postContent ) ) != "}" ) ) {
      print( paste( " php syntax error: ", postContent ) )
    } else {
      print( paste( " invalid JSON: ", postContent ) )
    }
    return( FALSE )
  }
  
  jsonData <- fromJSON( postContent )
  print( jsonData )
  
  if( "error" %in% names( jsonData ) ) {
    return( FALSE )
  }
  return( TRUE )
}

# example:
#   data <- wideDataSelect("lkpStreamCode", list("*"), orderBy=list("StreamCode DESC"))
#   exportCSV(data, "c:/temp/output.csv")

exportCSV<-function( data, filepath )
{
  if( getOption( "wideFunCall", 0 ) ) {
    print( call( "exportCSV", data = data, filepath = filepath ) )
  }
  
  # PARAMETERS (row.names, na)
  # - do not include row numbers
  # - save null values as conventional nulls instead of empty strings (i.e., "column1","column2",,"column4" instead of "column1","column2","","column4")
  write.csv( data, filepath, row.names = FALSE, na = "" )
}

# to read a csv, see read.csv(...)

1 Source

library(tidyverse)
library(sf)

# Columns that uniquely identify observations
ID_cols <- c(
  "ID",
  "StreamName",
  "StreamCode",
  "SiteCode",
  "Latitude",
  "Longitude",
  "SampleDate",
  "Sample",
  'OrgNameList'
)

# Endpoint names and groups
met_cols <- list(
  `Morphology` = list(
    `Width - Mean` = "AvgWettedWidth",
    `Depth - Mean` = "meanDepth",
    `Width:Depth Ratio` = "widthDepthRatio",
    `Usable Area` = "UsableArea"
  ),
  `Substrate Sizes` = list(
    `D15` = "Substd15",
    `D50` = "Substd50",
    `D85` = "Substd85",
    `D50Max` = "substPartMaxd50"
  ),
  `Indices` = list(
    #`Substrate Quality Rating` = "Substrate Quality Rating",
    `Sediment Sorting Index` = "sedimentSortingIndex",
    `Sediment Transport Potential` = "sedimentTransportPotential",
    `Bank Stability - Mean` = "bankStabilityAVG",
    `Habitat Homogeneity` = "scoreHomogeneity",
    `Habitat Stability Score` = "HabitatStabilityScore"
  ),
  `Percent` = list(
    `Filamentous Algae Percent` = "TotalVegTypePCFL",
    `Non-Filamentous Algae Percent` = "TotalVegTypePCAL",
    `Moss Percent` = "TotalVegTypePCSS",
    `Macrophytes Percent` = "TotalVegTypePCMC"
  ),
  `Cover` = list(
    `Wood Cover - Unembedded` = "UnembCoverTypePCWood",
    `Wood Cover - Total` = "TotalCoverTypePCWood",
    `Round Rock Cover - Unembedded` = "UnembCoverTypePCRndRock",
    `Round Rock Cover - Total` = "TotalCoverTypePCRndRock",
    `Flat Rock Cover - Unembedded` = "UnembCoverTypePCFltRock",
    `Flat Rock Cover - Total` = "TotalCoverTypePCFltRock",
    `Macrophyte Cover - Unembedded` = "UnembCoverTypePCMacrophyte",
    `Macrophyte Cover - Total` = "TotalCoverTypePCMacrophyte",
    `Bank Cover - Unembedded` = "UnembCoverTypePCBank",
    `Bank Cover - Total` = "TotalCoverTypePCBank"
  )
)

# Selected Fish Endpoints
sel_fish_ep<-c(
  `Community Total Biomass`="Comm_Biomass",
  `Community Total Density`="Comm_Abundance",
  `Community Total Richness`="Comm_Richness",
  `SATI Biomass`="Therm_SATI_Bioperc",
  `Cold Water Stenotherm Biomass`="Therm_cold_Bio",
  `% Cold Water Stenotherm Biomass`="Therm_cold_Bioperc",
  `Tolerant Taxa Biomass`="Tol_tol_Bio",
  `% Tolerant Taxa Biomass`="Tol_tol_Bioperc",
  `Intolerant Taxa Biomass`="Tol_intol_Bio",
  `% Intolerant Taxa Biomass`="Tol_intol_Bioperc",
  
  `Tolerant Taxa Richness`="Tol_tol_Rich",
  `Intolerant Taxa Richness`="Tol_intol_Rich",
  
  `Sensitive Darters Biomass`="Spec_DRT_Bio",
  `% Sensitive Darters Biomass`="Spec_DRT_Bioperc",
  `Brook Trout Biomass`="Spec_BKT_Bio",
  `% Brook Trout Biomass`="Spec_BKT_Bioperc",
  `Creek Chub Biomass`="Spec_CKC_Bio",
  `% Creek Chub Biomass`="Spec_CKC_Bioperc",
  `Open Substratum/Substratum Choosers Biomass`="Rep_OSS_Bio",
  `% Open Substratum/Substratum Choosers Biomass`="Rep_OSS_Bioperc",
  `Lithophilic Spawners Biomass`="Rep_Lith_Bio",
  `% Lithophilic Spawners Biomass`="Rep_Lith_Bioperc"
  
)

# function to reformat model terms

get_terms <- function(mod,add_fixed_var=T){
  fe_tbl <- as.data.frame(glmmTMB::fixef(mod)$cond)
  colnames(fe_tbl)<-"coef"
  fe_tbl$se <- diag(glmmTMB::sandwich(mod))
  
  add_var <- fe_tbl["(Intercept)","se"]
  if (!add_fixed_var) add_var <- 0
  
  re_tbl <- as.data.frame(glmmTMB::ranef(mod, condVar = T)) %>% 
    filter(component == "cond") %>% 
    mutate(val=condval,
           sd=condsd)
  
  comp_sel <- unique(re_tbl$grpvar)
  u_comp_sel <- str_split(comp_sel,"\\:")
  u_comp_sel <- unique(unlist(u_comp_sel))
  u_comp_sel <- unique(gsub("`","",u_comp_sel))
  
  for (i in u_comp_sel) {
    comp_sel_sub <- comp_sel[sapply(comp_sel,function(x) any(grepl(i,x)))]
    re_tbl[[i]]<-NA_character_
    
    for (ii in comp_sel_sub) {
      pos <- str_split(ii,"\\:",simplify = T)
      pos <- gsub("`","",pos)
      pos <- which(pos==i)
      
      re_tbl[re_tbl$grpvar==ii,i] <- str_split(re_tbl[re_tbl$grpvar==ii,"grp"],"\\:",simplify= T)[,pos]
    }
  }
  
  for (i in unique(re_tbl$SiteCode[!is.na(re_tbl$SiteCode)])){
    
    for (ii in unique(re_tbl$term)){
      val_sum <- re_tbl$val[re_tbl$SiteCode %in% i & re_tbl$term %in% ii]
      sd_sum <- re_tbl$sd[re_tbl$SiteCode %in% i & re_tbl$term %in% ii]
      
      for (iii in unique(re_tbl$grpvar)[!grepl("SiteCode",unique(re_tbl$grpvar))]){
        val_sum <- c(val_sum, re_tbl$condval[re_tbl$grpvar %in% iii &
                                               re_tbl$grp %in% gsub(paste0(":",i),"",unique(re_tbl[["grp"]][re_tbl$SiteCode %in% i])) &
                                               re_tbl$term %in% ii
        ]
        )
        sd_sum <- c(sd_sum, re_tbl$condsd[re_tbl$grpvar %in% iii &
                                            re_tbl$grp %in% gsub(paste0(":",i),"",unique(re_tbl[["grp"]][re_tbl$SiteCode %in% i])) &
                                            re_tbl$term %in% ii
        ]
        )
      }
      
      re_tbl$val[re_tbl$SiteCode %in% i & re_tbl$term %in% ii] <- sum(val_sum)
      re_tbl$sd[re_tbl$SiteCode %in% i & re_tbl$term %in% ii] <- sqrt(sum(sd_sum^2)+add_var^2)
      
    }
  }
  
  comp_list <- re_tbl %>% 
    filter(!is.na(SiteCode))%>% 
    mutate(group_var="SiteCode") %>% 
    select(group_var,group=SiteCode,component,term,val,sd,any_of(u_comp_sel)) %>% 
    group_by(group_var,group,component,term) %>% 
    summarise(
      value = sum(val),
      sd = sqrt(sum(sd^2)+add_var^2),
      across(any_of(u_comp_sel),~unique(.x)),
      .groups = "drop"
    )
  
  comp_list <- bind_rows(
    comp_list,
    re_tbl %>% 
      filter(is.na(SiteCode)) %>% 
      filter(!grpvar %in% unique(comp_list$group_var)) %>% 
      select(group_var=grpvar,
             group=grp,
             component,
             term,
             val=condval,
             sd=condsd,
             any_of(u_comp_sel)) %>% 
      group_by(across(c(any_of(u_comp_sel),group_var,group,component,term))) %>% 
      summarise(
        value = sum(val),
        sd = sqrt(sum(sd^2)+add_var^2),
        .groups = "drop"
      )
  )
  
  comp_list <- comp_list %>% 
    mutate(
      pred_low = value - (qnorm(1 - 0.1 / 2) * sd),
      pred_high = value + (qnorm(1 - 0.1 / 2) * sd),
      pred_low_trend = value - (qnorm(1 - 0.5 / 2) * sd),
      pred_high_trend = value + (qnorm(1 - 0.5 / 2) * sd),
    ) 
  
  return(comp_list)
  
}

get_terms2 <- function(mod,add_fixed_var=F,terms){

  fixed_list <- list(
    #"`UCA Class`:`Stress Class`:`Habitat Stability Class`" = as.formula(~`UCA Class`*`Stress Class`*`Habitat Stability Class`),
    #"`UCA Class`:`Habitat Stability Class`" = as.formula(~`UCA Class`*`Habitat Stability Class`),
    "`UCA Class`:`Stress Class`" = as.formula(~`UCA Class`*`Stress Class`),
    #"`Habitat Stability Class`" = as.formula(~`Habitat Stability Class`),
    "`Stress Class`" = as.formula(~`Stress Class`),
    "`UCA Class`" = as.formula(~`UCA Class`),
    "`1`" = as.formula(~1)
  )
  
  fixed_list <- map(fixed_list,
                    ~{
                      emt95 <- emmeans::emtrends(mod,.x,var="Year_n",level=0.95) %>% # consider: addl.vars = "SiteCode" ,component="response"
                        data.frame(check.names = F) %>% 
                        mutate(level=0.95,term="Year_n") %>% 
                        rename(val=Year_n.trend)
                      
                      emt50 <- emmeans::emtrends(mod,.x,var="Year_n",level=0.5) %>% 
                        data.frame(check.names = F) %>% 
                        mutate(level=0.5,term="Year_n") %>% 
                        rename(val=Year_n.trend)
                      
                      emm95 <- emmeans::emmeans(mod,.x,level=0.95, at = list(Year_n = max(mod$frame$Year_n)),type="response") %>% 
                        data.frame(check.names = F) %>% 
                        mutate(level=0.95,term="(Intercept)") %>% 
                        rename(val=response)
                      emm50 <- emmeans::emmeans(mod,.x,level=0.50, at = list(Year_n = max(mod$frame$Year_n)),type="response") %>% 
                        data.frame(check.names = F) %>% 
                        mutate(level=0.5,term="(Intercept)") %>% 
                        rename(val=response)
                      
                      bind_rows(emm95,emm50,emt95,emt50) %>% 
                        rowwise() %>% 
                        mutate(grp=paste0(c_across(any_of(c("UCA Class","Stress Class"))),collapse = ":"))
                    })
  
  fixed_list <- bind_rows(fixed_list,.id="grpvar") %>% 
    mutate(component="fixed") 
  
  re_tbl <- as.data.frame(glmmTMB::ranef(mod, condVar = T)) %>% 
    filter(component == "cond") %>% 
    filter(grpvar == "`UCA Class`:`Stress Class`:SiteCode") %>% 
    mutate(SiteCode = str_split(grp,":",simplify = T)[,3]) %>% 
    rename(val=condval,
           SE=condsd)
  
  out_frame_raw <- distinct(mod$frame[,sapply(mod$frame,class)=="factor"]) %>% 
    mutate(val = NA_real_,
           SE = NA_real_)
  
  out_frame <- out_frame_raw
  
  for (i in unique(out_frame$SiteCode)) {
    sub_frame1 <- fixed_list %>% 
      filter(term == "Year_n") %>% 
      filter(level == 0.95) %>% 
      filter(
        `UCA Class` == out_frame[out_frame$SiteCode==i,"UCA Class"],
        #`Habitat Stability Class` == out_frame[out_frame$SiteCode==i,"Habitat Stability Class"],
        `Stress Class` == out_frame[out_frame$SiteCode==i,"Stress Class"]
      )
    
    sub_frame2 <- re_tbl %>% 
      filter(SiteCode==i) %>% 
      filter(term == "Year_n")
    
    out_frame$val[out_frame$SiteCode==i] <- sum(sub_frame1$val,na.rm = T) + sum(sub_frame2$val,na.rm = T)
    out_frame$SE[out_frame$SiteCode==i] <- sqrt(sum(sub_frame1$SE^2,na.rm = T) + sum(sub_frame2$SE^2,na.rm = T))
    
  }
  
  out_frame <- out_frame %>% 
    mutate(
      pred_low = val - (qnorm(1 - 0.1 / 2) * SE),
      pred_high = val + (qnorm(1 - 0.1 / 2) * SE),
      pred_low_trend = val - (qnorm(1 - 0.5 / 2) * SE),
      pred_high_trend = val + (qnorm(1 - 0.5 / 2) * SE),
    ) %>% 
    mutate(
      term = "Year_n",
      component = "cond",
      grpvar = "`UCA Class`:`Stress Class`:SiteCode",
    )
  
  fixed_list_wide <- fixed_list %>% 
    pivot_wider(names_from = level, values_from = c(asymp.LCL,asymp.UCL)) %>% 
    rename(
      pred_low = asymp.LCL_0.95,
      pred_low_trend = asymp.LCL_0.5,
      pred_high = asymp.UCL_0.95,
      pred_high_trend = asymp.UCL_0.5
    )
  
  
  
  return(
    bind_rows(out_frame,fixed_list_wide) %>% 
      rename(value=val,
             group = grp,
             group_var = grpvar,
             sd = SE)
  )
  
}

summary_fun <- function(x){
  out <- map(
    c("Subshed","WSHED","UCA Class","Stress Class","UCA Class:Stress Class","UCA Class:Stress Class"),#,"UCA Class:Habitat Stability Class","UCA Class:Stress Class:Habitat Stability Class"
    ~{x %>% 
        as_tibble() %>% 
        mutate(`UCA Class:Stress Class` = paste(`UCA Class`,`Stress Class`,sep = ".")) %>% 
        #mutate(`UCA Class:Habitat Stability Class` = paste(`UCA Class`,`Habitat Stability Class`,sep = ".")) %>% 
        #mutate(`UCA Class:Stress Class:Habitat Stability Class` = paste(`UCA Class`,`Stress Class`,`Habitat Stability Class`,sep = ".")) %>% 
        mutate(group=.x) %>% 
        rename_with(.cols=any_of(.x),~paste0("group_val")) %>% 
        mutate(group_val=as.character(group_val)) %>% 
        group_by(group,group_val) %>% 
        summarise(
          across(c(
            `Upstream Catchment Area`,
          ),
          ~mean(log10(.x),na.rm=T)^10
          ),
          across(c(
            `Reach Slope`,
          ),
          ~expm1(mean(log1p(.x),na.rm=T))
          ),
          across(c(
            `Baseflow Index`,
            `Upstream catchment developed (%)`,
            `Upstream catchment agricultural (%)`,
            `Upstream catchment stress index`
          ),
          ~median(.x,na.rm=T)
          ),
          .groups="drop"
        )
    }
  )
  
  out <- bind_rows(out)
  return(out)
}

# Function to extract model optim

glmmTMB_get_optimum <- function(mod) {
  # Fixed effects
  starting_point <- list(
    beta = glmmTMB::fixef(mod)$cond,
    betazi = glmmTMB::fixef(mod)$zi,
    betadisp = glmmTMB::fixef(mod)$disp
  )
  #Random effects
  ref <- glmmTMB::ranef(mod)
  if (length(ref$cond) > 0) {
    starting_point$theta = ref$cond
  }
  if (length(ref$zi) > 0) {
    starting_point$thetazi = ref$zi
  }
  if (length(ref$disp) > 0) {
    starting_point$thetadisp = ref$disp
  }
  return(starting_point)
}

# Function to generate heatmaps
heatmap_fn<-function(x_df,
                     y_df=NULL,
                     alpha=0.1,
                     sig_est=0.6,
                     method="spearman",
                     scale_grad=ggplot2::scale_fill_gradient2(
                       low = scales::muted(rgb(red=0/255,green=102/255,blue=255/255)),
                       mid = "white",
                       high = scales::muted(rgb(red=0/255,green=176/255,blue=80/255)),
                       limits=c(-1,1)
                     ),
                     ...){
  
  if (is.null(y_df)){
    y_df<-x_df
  }
  cor_list<-lapply(x_df,function(x) lapply(y_df,function(y) broom::tidy(cor.test(x,y,method=method,...))))
  
  pcorr_out<-purrr::map_dfr(cor_list,.id="x_var",~purrr::map_dfr(.,dplyr::bind_rows,.id="y_var"))
  
  pcorr_out <- pcorr_out %>%
    mutate(p.value = case_when(
      estimate >= 1 - .Machine$double.eps ~ NA,
      estimate <= -1 + .Machine$double.eps~ NA,
      T ~ p.value
    )) %>%
    mutate(estimate = case_when(
      estimate >= 1 - .Machine$double.eps ~ NA,
      estimate <= -1 + .Machine$double.eps~ NA,
      T ~ estimate
    )) 
  
  rho<-pcorr_out %>%
    select(x_var,y_var,estimate) %>%
    spread(x_var,estimate) %>%
    data.frame(check.names = F) %>%
    tibble::column_to_rownames("y_var") %>%
    textshape::cluster_matrix(method ="ward.D2")
  
  pv<-pcorr_out %>%
    mutate(p.value=case_when(
      is.na(p.value) ~ "darkgray",
      p.value>alpha ~ "darkgray",
      p.value<=alpha & abs(estimate)>=sig_est ~ "darkred",
      p.value<=alpha ~ "orange",
      T ~ "darkgray"
    )) %>%
    select(x_var,y_var,p.value) %>%
    spread(x_var,p.value) %>%
    data.frame(check.names = F) %>%
    tibble::column_to_rownames("y_var")
  
  pv<-pv[rownames(rho),colnames(rho)]
  
  plot_out <- pcorr_out %>%
    select(x_var,y_var,estimate) %>% 
    spread(x_var,estimate) %>%
    data.frame(check.names = F) %>%
    tibble::column_to_rownames("y_var") %>%
    heatmaply::ggheatmap(hclust_method ="ward.D2",
                         seriate="GW",
                         margins =c(40,40,40,40),
                         branches_lwd = 0.25,
                         grid_gap=3,
                         heatmap_layers=list(
                           theme(
                             axis.text.y = element_text(size = 8)   ,
                             axis.text.x = element_text(
                               size = 8,
                               hjust = 1,
                               vjust = 1,
                               angle = 45
                             ),
                             legend.position = "none"),
                           geom_tile(data=pv %>%
                                       rownames_to_column("Endpoint") %>%
                                       gather("Parameter","Value",-Endpoint),
                                     aes(x=Parameter,y=Endpoint),
                                     #alpha=0,
                                     colour=pv %>%
                                       rownames_to_column("Endpoint") %>%
                                       gather("Parameter","Value",-Endpoint) %>%
                                       pull(Value),
                                     width=0.975, height=0.975,
                                     size=0.3,
                                     fill=NA,
                                     inherit.aes = F),
                           guides(colour="none")
                         ),
                         scale_fill_gradient_fun = scale_grad)
  
  cowplot::ggdraw(plot_out)
  
}

# transformer for car::logit
clogit <- scales::new_transform(
  name="clogit",
  transform = function(x) {
    x[x<0] <- 0
    x[x>1] <- 1
    car::logit(x,T,0.025)
  },
  inverse = function(x) {
    a <- 0.025
    a <- (1-2*a)
    a <- (a*(1+exp(x))+(exp(x)-1))/(2*a*(1+exp(x)))
    zapsmall(a)
  },
  domain = c(0, 1)
)

# Functions to plot axis on pseudo-log scales
force_all <- function(...) list(...) # Helper function defined in scales

pseudo_log <- function (val, sigma = 1, base = exp(1)) {
  asinh(val/(2*sigma))/log(base)
}

pseudo_log_breaks <- function(n = 5, sigma = 1, base = 10) {
  force_all(n, base)
  n_default = n
  
  force_all <- function(...) list(...) # Helper function defined in scales
  
  pseudo_log <- function (val, sigma = 1, base = exp(1)) {
    asinh(val/(2*sigma))/log(base)
  }
  
  function(x, n = n_default) {
    rng <- pseudo_log(range(x, na.rm = TRUE), sigma = sigma, base = base)
    min <- floor(rng[1])
    max <- ceiling(rng[2])
    
    if (max == min) return(base^min)
    
    by <- floor((max - min) / n) + 1
    breaks <- base^seq(min, max, by = by)
    relevant_breaks <- base^rng[1] <= breaks & breaks <= base^rng[2]
    if (sum(relevant_breaks) >= (n - 2)) return(breaks)
    
    # the easy solution to get more breaks is to decrease 'by'
    while (by > 1) {
      by <- by - 1
      breaks <- base^seq(min, max, by = by)
      relevant_breaks <- base^rng[1] <= breaks & breaks <= base^rng[2]
      if (sum(relevant_breaks) >= (n - 2)) return(breaks)
    }
    scales:::log_sub_breaks(rng, n = n, base = base)
  }
}

# Helper function to extract a legend
get_legend <- function(p, position = NULL){
  
  if(.is_list(p)){
    continue <- TRUE
    i <- 1
    while(i <= length(p) & continue){
      leg <- .get_legend(p[[i]], position = position)
      if(!is.null(leg)) continue <- FALSE
      i <- i+1
    }
  }
  else{
    leg <- .get_legend(p, position = position)
  }
  leg
}

.is_list <- function(x){
  inherits(x, "list")
}


# Return legend for one plot
.get_legend <- function(p, position = NULL){
  
  if(is.null(p)) return(NULL)
  if(!is.null(position)){
    p <- p + theme(legend.position = position)
  }
  tmp <- ggplot_gtable(ggplot_build(p))
  leg <- which(sapply(tmp$grobs, function(x) x$name) == "guide-box")
  if(length(leg) > 0) leg <- tmp$grobs[[leg]]
  else leg <- NULL
  leg
}

col2hex <- function(color, alpha = 1) {
  
  if (alpha < 0 | alpha > 1)
    stop('alpha must be between 0 and 1.')
  
  if (missing(color))
    stop('Need a color')
  
  if (!is.character(color))
    stop(strwrap('Elements must be character string as a named R color or
         hex (e.g. "#ffffff")'))
  
  if (!color %in% colors()) {
    if (!nchar(color) == 7 | !grepl('^#', color))
      stop(strwrap('color must be an R color, e.g. see colors(),
                   or a hex of the form #ff5500'))
  }
  
  rgb(t(col2rgb(color)/255), alpha = alpha)
}

# Calculate fish endpoints
calc_fish_ep<-function(df) {
  
  ep_tbl<-df %>% 
    group_by(SampleEventID,SampleDate) %>% 
    mutate(
      # General Community
      Comm_Biomass=sum(EstimatedBiomass,na.rm=T),
      Comm_Abundance=sum(NumberOfFish,na.rm=T),
      Comm_Biomass_SATI=sum(EstimatedBiomass[!is.na(OPT_TEMP)],na.rm=T),
      Comm_Abundance_SATI=sum(NumberOfFish[!is.na(OPT_TEMP)],na.rm=T),
      Comm_Richness=length(unique(SpeciesCode)),
      EstimatedBiomass_perc=EstimatedBiomass/Comm_Biomass,
      NumberOfFish_perc=NumberOfFish/Comm_Abundance,
      EstimatedBiomass_perc_SATI=EstimatedBiomass/Comm_Biomass_SATI,
      NumberOfFish_perc_SATI=NumberOfFish/Comm_Abundance_SATI,
      SATI_EstimatedBiomass_perc=EstimatedBiomass_perc_SATI*OPT_TEMP,
      SATI_NumberOfFish_perc=NumberOfFish_perc_SATI*OPT_TEMP,
    ) %>% 
    summarise(
      # General Community
      Comm_Biomass=head(Comm_Biomass,1),
      Comm_Abundance=head(Comm_Abundance,1),
      Comm_Richness=head(Comm_Richness,1),
      # Thermal Regime 
      Therm_SATI_Bioperc=sum(SATI_EstimatedBiomass_perc,na.rm=T),
      Therm_SATI_Abuperc=sum(SATI_NumberOfFish_perc,na.rm=T),
      
      Therm_cold_Bioperc=sum(EstimatedBiomass_perc[Thermal.Regime=="coldwater"],na.rm=T),
      Therm_cool_Bioperc=sum(EstimatedBiomass_perc[Thermal.Regime=="coolwater"],na.rm=T),
      Therm_coolcold_Bioperc=sum(EstimatedBiomass_perc[Thermal.Regime %in% c("coolwater","coldwater")],na.rm=T),
      Therm_warm_Bioperc=sum(EstimatedBiomass_perc[Thermal.Regime %in% c("warmwater")],na.rm=T),
      Therm_cold_Abuperc=sum(NumberOfFish_perc[Thermal.Regime=="coldwater"],na.rm=T),
      Therm_cool_Abuperc=sum(NumberOfFish_perc[Thermal.Regime=="coolwater"],na.rm=T),
      Therm_coolcold_Abuperc=sum(NumberOfFish_perc[Thermal.Regime %in% c("coolwater","coldwater")],na.rm=T),
      Therm_warm_Abuperc=sum(NumberOfFish_perc[Thermal.Regime %in% c("warmwater")],na.rm=T),
      Therm_cold_Bio=sum(EstimatedBiomass_perc[Thermal.Regime=="coldwater"],na.rm=T)*Comm_Biomass,
      Therm_cool_Bio=sum(EstimatedBiomass_perc[Thermal.Regime=="coolwater"],na.rm=T)*Comm_Biomass,
      Therm_coolcold_Bio=sum(EstimatedBiomass_perc[Thermal.Regime %in% c("coolwater","coldwater")],na.rm=T)*Comm_Biomass,
      Therm_warm_Bio=sum(EstimatedBiomass_perc[Thermal.Regime %in% c("warmwater")],na.rm=T)*Comm_Biomass,
      Therm_cold_Abu=sum(NumberOfFish_perc[Thermal.Regime=="coldwater"],na.rm=T)*Comm_Abundance,
      Therm_cool_Abu=sum(NumberOfFish_perc[Thermal.Regime=="coolwater"],na.rm=T)*Comm_Abundance,
      Therm_coolcold_Abu=sum(NumberOfFish_perc[Thermal.Regime %in% c("coolwater","coldwater")],na.rm=T)*Comm_Abundance,
      Therm_warm_Abu=sum(NumberOfFish_perc[Thermal.Regime %in% c("warmwater")],na.rm=T)*Comm_Abundance,
      Therm_warm_Rich=length(NumberOfFish_perc[Thermal.Regime %in% c("warmwater")]),
      Therm_cold_Rich=length(NumberOfFish_perc[Thermal.Regime %in% c("coldwater")]),
      Therm_cool_Rich=length(NumberOfFish_perc[Thermal.Regime %in% c("coolwater")]),
      Therm_coolcold_Rich=length(NumberOfFish_perc[Thermal.Regime %in% c("coolwater","coldwater")]),
      # Trophic Guild 
      Troph_invert_Bioperc=sum(EstimatedBiomass_perc[grepl("invertivore",Trophic.Class)],na.rm=T),
      Troph_detrit_Bioperc=sum(EstimatedBiomass_perc[grepl("detritivore",Trophic.Class)],na.rm=T),
      Troph_carn_Bioperc=sum(EstimatedBiomass_perc[grepl("carnivore",Trophic.Class)],na.rm=T),
      Troph_herb_Bioperc=sum(EstimatedBiomass_perc[grepl("herbivore",Trophic.Class)],na.rm=T),
      Troph_plank_Bioperc=sum(EstimatedBiomass_perc[grepl("planktivore",Trophic.Class)],na.rm=T),
      Troph_invert_Bio=sum(EstimatedBiomass_perc[grepl("invertivore",Trophic.Class)],na.rm=T)*Comm_Biomass,
      Troph_detrit_Bio=sum(EstimatedBiomass_perc[grepl("detritivore",Trophic.Class)],na.rm=T)*Comm_Biomass,
      Troph_carn_Bio=sum(EstimatedBiomass_perc[grepl("carnivore",Trophic.Class)],na.rm=T)*Comm_Biomass,
      Troph_herb_Bio=sum(EstimatedBiomass_perc[grepl("herbivore",Trophic.Class)],na.rm=T)*Comm_Biomass,
      Troph_plank_Bio=sum(EstimatedBiomass_perc[grepl("planktivore",Trophic.Class)],na.rm=T)*Comm_Biomass,
      Troph_invert_Abuperc=sum(NumberOfFish_perc[grepl("invertivore",Trophic.Class)],na.rm=T),
      Troph_detrit_Abuperc=sum(NumberOfFish_perc[grepl("detritivore",Trophic.Class)],na.rm=T),
      Troph_carn_Abuperc=sum(NumberOfFish_perc[grepl("carnivore",Trophic.Class)],na.rm=T),
      Troph_herb_Abuperc=sum(NumberOfFish_perc[grepl("herbivore",Trophic.Class)],na.rm=T),
      Troph_plank_Abuperc=sum(NumberOfFish_perc[grepl("planktivore",Trophic.Class)],na.rm=T),
      Troph_invert_Abu=sum(NumberOfFish_perc[grepl("invertivore",Trophic.Class)],na.rm=T)*Comm_Abundance,
      Troph_detrit_Abu=sum(NumberOfFish_perc[grepl("detritivore",Trophic.Class)],na.rm=T)*Comm_Abundance,
      Troph_carn_Abu=sum(NumberOfFish_perc[grepl("carnivore",Trophic.Class)],na.rm=T)*Comm_Abundance,
      Troph_herb_Abu=sum(NumberOfFish_perc[grepl("herbivore",Trophic.Class)],na.rm=T)*Comm_Abundance,
      Troph_plank_Abu=sum(NumberOfFish_perc[grepl("planktivore",Trophic.Class)],na.rm=T)*Comm_Abundance,
      Troph_invert_Rich=length(NumberOfFish_perc[grepl("invertivore",Trophic.Class)]),
      Troph_detrit_Rich=length(NumberOfFish_perc[grepl("detritivore",Trophic.Class)]),
      Troph_carn_Rich=length(NumberOfFish_perc[grepl("carnivore",Trophic.Class)]),
      Troph_herb_Rich=length(NumberOfFish_perc[grepl("herbivore",Trophic.Class)]),
      Troph_plank_Rich=length(NumberOfFish_perc[grepl("planktivore",Trophic.Class)]),
      # Tolerance
      Tol_tol_Bioperc=sum(EstimatedBiomass_perc[Tolerance=="tolerant"],na.rm=T),
      Tol_intol_Bioperc=sum(EstimatedBiomass_perc[Tolerance=="intolerant"],na.rm=T),
      Tol_tol_Abuperc=sum(NumberOfFish_perc[Tolerance=="tolerant"],na.rm=T),
      Tol_intol_Abuperc=sum(NumberOfFish_perc[Tolerance=="intolerant"],na.rm=T),
      Tol_tol_Bio=sum(EstimatedBiomass_perc[Tolerance=="tolerant"],na.rm=T)*Comm_Biomass,
      Tol_intol_Bio=sum(EstimatedBiomass_perc[Tolerance=="intolerant"],na.rm=T)*Comm_Biomass,
      Tol_tol_Abu=sum(NumberOfFish_perc[Tolerance=="tolerant"],na.rm=T)*Comm_Abundance,
      Tol_intol_Abu=sum(NumberOfFish_perc[Tolerance=="intolerant"],na.rm=T)*Comm_Abundance,
      Tol_tol_Rich=length(NumberOfFish_perc[Tolerance=="tolerant"]),
      Tol_intol_Rich=length(NumberOfFish_perc[Tolerance=="intolerant"]),
      # Individual Species - 
      Spec_RBD_Bioperc=sum(EstimatedBiomass_perc[SpeciesCode==337],na.rm=T),
      Spec_RBD_Abuperc=sum(NumberOfFish_perc[SpeciesCode==337],na.rm=T),
      Spec_BKT_Bioperc=sum(EstimatedBiomass_perc[SpeciesCode==80],na.rm=T),
      Spec_BKT_Abuperc=sum(NumberOfFish_perc[SpeciesCode==80],na.rm=T),
      Spec_GDF_Bioperc=sum(EstimatedBiomass_perc[SpeciesCode==181],na.rm=T),
      Spec_GDF_Abuperc=sum(NumberOfFish_perc[SpeciesCode==181],na.rm=T),
      Spec_CMC_Bioperc=sum(EstimatedBiomass_perc[SpeciesCode==186],na.rm=T),
      Spec_CMC_Abuperc=sum(NumberOfFish_perc[SpeciesCode==186],na.rm=T),
      Spec_RDG_Bioperc=sum(EstimatedBiomass_perc[SpeciesCode==366],na.rm=T),
      Spec_RDG_Abuperc=sum(NumberOfFish_perc[SpeciesCode==366],na.rm=T),
      Spec_CKC_Bioperc=sum(EstimatedBiomass_perc[SpeciesCode==212],na.rm=T),
      Spec_CKC_Abuperc=sum(NumberOfFish_perc[SpeciesCode==212],na.rm=T),
      Spec_SKP_Bioperc=sum(EstimatedBiomass_perc[SpeciesCode %in% c(381,382)],na.rm=T),
      Spec_SKP_Abuperc=sum(NumberOfFish_perc[SpeciesCode %in% c(381,382)],na.rm=T),
      Spec_LMP_Bioperc=sum(EstimatedBiomass_perc[SpeciesCode %in% c(11:14)],na.rm=T),
      Spec_LMP_Abuperc=sum(NumberOfFish_perc[SpeciesCode %in% c(11:14)],na.rm=T),
      Spec_DRT_Bioperc=sum(EstimatedBiomass_perc[SpeciesCode %in% c(335:340,342)],na.rm=T),
      Spec_DRT_Abuperc=sum(NumberOfFish_perc[SpeciesCode %in% c(335:340,342)],na.rm=T),
      Spec_CEB_Bioperc=sum(EstimatedBiomass_perc[SpeciesCode %in% c(311, 313:317)],na.rm=T),
      Spec_CEB_Abuperc=sum(NumberOfFish_perc[SpeciesCode %in% c(311, 313:317)],na.rm=T),
      Spec_CTF_Bioperc=sum(EstimatedBiomass_perc[SpeciesCode %in% c(231:233)],na.rm=T),
      Spec_CTF_Abuperc=sum(NumberOfFish_perc[SpeciesCode %in% c(231:233)],na.rm=T),
      Spec_SAL_Bioperc=sum(EstimatedBiomass_perc[SpeciesCode %in% c(71:83)],na.rm=T),
      Spec_SAL_Abuperc=sum(NumberOfFish_perc[SpeciesCode %in% c(71:83)],na.rm=T),
      
      Spec_RBD_Bio=sum(EstimatedBiomass_perc[SpeciesCode==337],na.rm=T)*Comm_Biomass,
      Spec_RBD_Abu=sum(NumberOfFish_perc[SpeciesCode==337],na.rm=T)*Comm_Abundance,
      Spec_BKT_Bio=sum(EstimatedBiomass_perc[SpeciesCode==80],na.rm=T)*Comm_Biomass,
      Spec_BKT_Abu=sum(NumberOfFish_perc[SpeciesCode==80],na.rm=T)*Comm_Abundance,
      Spec_GDF_Bio=sum(EstimatedBiomass_perc[SpeciesCode==181],na.rm=T)*Comm_Biomass,
      Spec_GDF_Abu=sum(NumberOfFish_perc[SpeciesCode==181],na.rm=T)*Comm_Abundance,
      Spec_CMC_Bio=sum(EstimatedBiomass_perc[SpeciesCode==186],na.rm=T)*Comm_Biomass,
      Spec_CMC_Abu=sum(NumberOfFish_perc[SpeciesCode==186],na.rm=T)*Comm_Abundance,
      Spec_RDG_Bio=sum(EstimatedBiomass_perc[SpeciesCode==366],na.rm=T)*Comm_Biomass,
      Spec_RDG_Abu=sum(NumberOfFish_perc[SpeciesCode==366],na.rm=T)*Comm_Abundance,
      Spec_CKC_Bio=sum(EstimatedBiomass_perc[SpeciesCode==212],na.rm=T)*Comm_Biomass,
      Spec_CKC_Abu=sum(NumberOfFish_perc[SpeciesCode==212],na.rm=T)*Comm_Abundance,
      Spec_SKP_Bio=sum(EstimatedBiomass_perc[SpeciesCode %in% c(381,382)],na.rm=T)*Comm_Biomass,
      Spec_SKP_Abu=sum(NumberOfFish_perc[SpeciesCode %in% c(381,382)],na.rm=T)*Comm_Abundance,
      Spec_LMP_Bio=sum(EstimatedBiomass_perc[SpeciesCode %in% c(11:14)],na.rm=T)*Comm_Biomass,
      Spec_LMP_Abu=sum(NumberOfFish_perc[SpeciesCode %in% c(11:14)],na.rm=T)*Comm_Abundance,
      Spec_DRT_Bio=sum(EstimatedBiomass_perc[SpeciesCode %in% c(335:340,342)],na.rm=T)*Comm_Biomass,
      Spec_DRT_Abu=sum(NumberOfFish_perc[SpeciesCode %in% c(335:340,342)],na.rm=T)*Comm_Abundance,
      Spec_CEB_Bio=sum(EstimatedBiomass_perc[SpeciesCode %in% c(311, 313:317)],na.rm=T)*Comm_Biomass,
      Spec_CEB_Abu=sum(NumberOfFish_perc[SpeciesCode %in% c(311, 313:317)],na.rm=T)*Comm_Abundance,
      Spec_CTF_Bio=sum(EstimatedBiomass_perc[SpeciesCode %in% c(231:233)],na.rm=T)*Comm_Biomass,
      Spec_CTF_Abu=sum(NumberOfFish_perc[SpeciesCode %in% c(231:233)],na.rm=T)*Comm_Abundance,
      Spec_SAL_Bio=sum(EstimatedBiomass_perc[SpeciesCode %in% c(71:83)],na.rm=T)*Comm_Biomass,
      Spec_SAL_Abu=sum(NumberOfFish_perc[SpeciesCode %in% c(71:83)],na.rm=T)*Comm_Abundance,
      # Environment
      Env_benthic_Bioperc=sum(EstimatedBiomass_perc[Environment=="benthic"],na.rm=T),
      Env_pelagic_Bioperc=sum(EstimatedBiomass_perc[Environment=="pelagic"],na.rm=T),
      Env_benthopelagic_Bioperc=sum(EstimatedBiomass_perc[Environment=="benthopelagic"],na.rm=T),
      Env_benthic_Abuperc=sum(NumberOfFish_perc[Environment=="benthic"],na.rm=T),
      Env_pelagic_Abuperc=sum(NumberOfFish_perc[Environment=="pelagic"],na.rm=T),
      Env_benthopelagic_Abuperc=sum(NumberOfFish_perc[Environment=="benthopelagic"],na.rm=T),
      Env_benthic_Bio=sum(EstimatedBiomass_perc[Environment=="benthic"],na.rm=T)*Comm_Biomass,
      Env_pelagic_Bio=sum(EstimatedBiomass_perc[Environment=="pelagic"],na.rm=T)*Comm_Biomass,
      Env_benthopelagic_Bio=sum(EstimatedBiomass_perc[Environment=="benthopelagic"],na.rm=T)*Comm_Biomass,
      Env_benthic_Abu=sum(NumberOfFish_perc[Environment=="benthic"],na.rm=T)*Comm_Abundance,
      Env_pelagic_Abu=sum(NumberOfFish_perc[Environment=="pelagic"],na.rm=T)*Comm_Abundance,
      Env_benthopelagic_Abu=sum(NumberOfFish_perc[Environment=="benthopelagic"],na.rm=T)*Comm_Abundance,
      Env_benthic_Rich=length(NumberOfFish_perc[Environment=="benthic"]),
      Env_pelagic_Rich=length(NumberOfFish_perc[Environment=="pelagic"]),
      Env_benthopelagic_Rich=length(NumberOfFish_perc[Environment=="benthopelagic"]),
      # Reproductive Guild
      ## Guarders
      Rep_guard_Bioperc=sum(EstimatedBiomass_perc[grepl("^B\\.",Reproductive.Code)],na.rm=T),
      Rep_guard_Abuperc=sum(NumberOfFish_perc[grepl("^B\\.",Reproductive.Code)],na.rm=T),
      Rep_guard_Bio=sum(EstimatedBiomass_perc[grepl("^B\\.",Reproductive.Code)],na.rm=T)*Comm_Biomass,
      Rep_guard_Abu=sum(NumberOfFish_perc[grepl("^B\\.",Reproductive.Code)],na.rm=T)*Comm_Abundance,
      Rep_guard_Rich=length(NumberOfFish_perc[grepl("^B\\.",Reproductive.Code)]),
      ## Nesting?
      Rep_OSS_Bioperc=sum(EstimatedBiomass_perc[grepl("\\.1\\.",Reproductive.Code)],na.rm=T),
      Rep_OSS_Abuperc=sum(NumberOfFish_perc[grepl("\\.1\\.",Reproductive.Code)],na.rm=T),
      Rep_OSS_Bio=sum(EstimatedBiomass_perc[grepl("\\.1\\.",Reproductive.Code)],na.rm=T)*Comm_Biomass,
      Rep_OSS_Abu=sum(NumberOfFish_perc[grepl("\\.1\\.",Reproductive.Code)],na.rm=T)*Comm_Abundance,
      Rep_OSS_Rich=length(NumberOfFish_perc[grepl("\\.1\\.",Reproductive.Code)]),
      Rep_BrHid_Bioperc=sum(EstimatedBiomass_perc[grepl("\\.2\\.",Reproductive.Code)],na.rm=T),
      Rep_BrHid_Abuperc=sum(NumberOfFish_perc[grepl("\\.2\\.",Reproductive.Code)],na.rm=T),
      Rep_BrHid_Bio=sum(EstimatedBiomass_perc[grepl("\\.2\\.",Reproductive.Code)],na.rm=T)*Comm_Biomass,
      Rep_BrHid_Abu=sum(NumberOfFish_perc[grepl("\\.2\\.",Reproductive.Code)],na.rm=T)*Comm_Abundance,
      Rep_BrHid_Rich=length(NumberOfFish_perc[grepl("\\.2\\.",Reproductive.Code)]),
      ## Lithophilic
      Rep_Lith_Bioperc=sum(EstimatedBiomass_perc[grepl("Litho",Reproductive.Guild)],na.rm=T),
      Rep_Lith_Abuperc=sum(NumberOfFish_perc[grepl("Litho",Reproductive.Guild)],na.rm=T),
      Rep_Lith_Bio=sum(EstimatedBiomass_perc[grepl("Litho",Reproductive.Guild)],na.rm=T)*Comm_Biomass,
      Rep_Lith_Abu=sum(NumberOfFish_perc[grepl("Litho",Reproductive.Guild)],na.rm=T)*Comm_Abundance,
      Rep_Lith_Rich=length(NumberOfFish_perc[grepl("Litho",Reproductive.Guild)]),
      .groups="drop"
    ) 
  
  return(ep_tbl)
}

2 Download Data

source("00_Functions/FWIS_Functions.R")
library(tidyverse)

username<-Sys.getenv("FWIS_username")
password<-Sys.getenv("FWIS_password")

fwis_login<-wideLogin("http://www.comap.ca/fwis/wideR.php",username,password)


# Habitat Data ------------------------------------------------------------

hab_tbls<-c("tblChannelStructureSummary")

columns <- wideDataSelect(1,
                          "tblDataAPI",
                          list("tbldColName", "tbldActionsRdRegUser", "tbldPK", "tbldNulls", "tbldType", "tbldGSConstraint","tbldFilePath","tbldFileFolder"),
                          where=list(tbldUserTblName=c("=",hab_tbls)),
                          orderBy=list("tbldOrdering"))

data <- wideDataSelect( 1, table=hab_tbls, colsSelected=as.list(columns$tbldColName) )

write_csv(data,"data/tblChannelStructureSummary.csv")

# Fish Data ---------------------------------------------------------------

fish_tbls<-c("tblFishSummaryOfTotalCatches")

columns <- wideDataSelect(1,
                          "tblDataAPI",
                          list("tbldColName", "tbldActionsRdRegUser", "tbldPK", "tbldNulls", "tbldType", "tbldGSConstraint","tbldFilePath","tbldFileFolder"),
                          where=list(tbldUserTblName=c("=",fish_tbls)),
                          orderBy=list("tbldOrdering"))

data <- wideDataSelect( 1, table=fish_tbls, colsSelected=as.list(columns$tbldColName) )

write_csv(data,"data/tblFishSummaryOfTotalCatches.csv")

# Fish Habitat Data -------------------------------------------------------

fish_tbls<-c("lkpSpeciesHabitatRequirement")

columns <- wideDataSelect(1,
                          "tblDataAPI",
                          list("tbldColName", "tbldActionsRdRegUser", "tbldPK", "tbldNulls", "tbldType", "tbldGSConstraint","tbldFilePath","tbldFileFolder"),
                          where=list(tbldUserTblName=c("=",fish_tbls)),
                          orderBy=list("tbldOrdering"))

data <- wideDataSelect( 1, table=fish_tbls, colsSelected=as.list(columns$tbldColName) )

write_csv(data,"data/lkpSpeciesHabitatRequirement.csv")

# Benthic Invertebrate Data -----------------------------------------------

# Benthic Invertebrate Download -------------------------------------------

tables <- unique(wideDataSelect(1,"tblDataAPI", list("tbldUserTblName"), orderBy = list("tbldUserTblName")))

bic_tbls<-tables$tbldUserTblName[grepl("benthic",tolower(tables$tbldUserTblName)) & grepl("tbl",tolower(tables$tbldUserTblName))]

bic_all<-lapply(setNames(bic_tbls,bic_tbls),function(x){
  columns <- wideDataSelect(1,
                            "tblDataAPI",
                            list("tbldColName", "tbldActionsRdRegUser", "tbldPK", "tbldNulls", "tbldType", "tbldGSConstraint","tbldFilePath","tbldFileFolder"),
                            where=list(tbldUserTblName=c("=",x)),
                            orderBy=list("tbldOrdering"))
  
  wideDataSelect( 1, table=x, colsSelected=as.list(columns$tbldColName) )
  
})

saveRDS(bic_all,"data/BIC_ALL.rds")

obbn_tbls<-tables$tbldUserTblName[grepl("obbn",tolower(tables$tbldUserTblName)) & grepl("tbl",tolower(tables$tbldUserTblName))]

obbn_all<-lapply(setNames(obbn_tbls,obbn_tbls),function(x){
  columns <- wideDataSelect(1,
                            "tblDataAPI",
                            list("tbldColName", "tbldActionsRdRegUser", "tbldPK", "tbldNulls", "tbldType", "tbldGSConstraint","tbldFilePath","tbldFileFolder"),
                            where=list(tbldUserTblName=c("=",x)),
                            orderBy=list("tbldOrdering"))
  
  columns$tbldColName <- gsub("Â", "", columns$tbldColName)
  
  wideDataSelect( 1, table=x, colsSelected=as.list(columns$tbldColName) )
  
})

saveRDS(obbn_all,"data/OBBN_ALL.rds")

3 Format Raw Data

source("00_source.R")
library(openxlsx)


# Read in Data ------------------------------------------------------------
dat <- read_csv("data/tblChannelStructureSummary.csv", locale = locale(encoding = "latin1"))


# Add additional endpoints ------------------------------------------------
dat <- dat %>% 
  mutate(
    TotalCoverTypePCWood = EmbCoverTypePCWood + UnembCoverTypePCWood,
    TotalCoverTypePCRndRock = EmbCoverTypePCRndRock + UnembCoverTypePCRndRock,
    TotalCoverTypePCFltRock = EmbCoverTypePCFltRock + UnembCoverTypePCFltRock,
    TotalCoverTypePCMacrophyte = EmbCoverTypePCMacrophyte + UnembCoverTypePCMacrophyte,
    TotalCoverTypePCBank = EmbCoverTypePCBank + UnembCoverTypePCBank,
  ) %>% 
  mutate(
    TotalVegTypePCFL = EmbVegTypePCFL + UnembVegTypePCFL,
    TotalVegTypePCAL = EmbVegTypePCAL + UnembVegTypePCAL,
    TotalVegTypePCSS = EmbVegTypePCSS + UnembVegTypePCSS,
    TotalVegTypePCMC = EmbVegTypePCMC + UnembVegTypePCMC,
  )

# Exclude test data
dat <- dat %>% 
  select(all_of(c(unlist(ID_cols),unlist(met_cols),StreamType = "substrateQualityRatingCode"))) %>% 
  filter(!StreamCode %in% c("TEST")) %>% 
  filter(!grepl("Training",OrgNameList)) %>% 
  filter(SiteCode != "NCD21-1C") #%>% # TRCA site not within bounds of TRCA shape file
#filter(SiteCode != "PT003WM") %>%  # incorrect coordinates
#filter(SiteCode != "PT002WM") 
# filter(month(SampleDate) != "1",
#        month(SampleDate) != "2",
#        month(SampleDate) != "12")

# Remove samples without location information
dat <- dat  %>% 
  rename_with(~gsub("^.*?\\.","",.x)) %>% 
  mutate(Longitude = as.numeric(Longitude),
         Latitude = as.numeric(Latitude)) %>% 
  group_by(SiteCode) %>% 
  mutate(Longitude = median(Longitude,na.rm = T),
         Latitude = median(Latitude,na.rm = T),
  ) %>% 
  filter(!is.na(Longitude),
         !is.na(Latitude)) %>% 
  ungroup()

# Add utility columns
dat <- dat %>% 
  group_by(!!!syms(ID_cols[!ID_cols %in% c("ID","SampleDate","Sample","Longitude","Latitude")]))  %>% 
  mutate(
    Month = month(SampleDate),
    Year = year(SampleDate),
    `Sample Season` = case_when(
      Month %in% 3:5 ~ "Spring",
      Month %in% 6:8 ~ "Summer",
      Month %in% 9:11 ~ "Fall",
      T ~ "Winter"
    ),
    `Sample Years` = paste(sort(unique(Year)),collapse = ", "),
    `Sample Months` = paste(sort(unique(Month)),collapse = ", "),
    `Years of Data` = length(unique(Year)),
    `Year-Month`=format(SampleDate,"%Y-%m"),
    `Sample Year-Months` = paste(sort(unique(`Year-Month`)),collapse = ", "),
    `Sample Year-Season` = paste(sort(unique(`Year-Month`)),collapse = ", "),
    `TRCA` = OrgNameList == "Toronto and Region Conservation Authority",
    `Sample Season`=factor(`Sample Season`,levels=c("Sping","Summer","Fall","Winter"))
  ) %>% 
  filter(Year < 2025) %>% 
  ungroup() %>% 
  group_by(!!!syms(c(ID_cols[!ID_cols %in% c("ID","SampleDate","Sample","Longitude","Latitude")],"Year"))) %>% 
  mutate( 
    `Sampled Days within Year` = length(unique(SampleDate)),
    `Sampled Seasons within Year` = length(unique(`Sample Season`)),
  ) %>% 
  ungroup() %>% 
  rowwise() %>% 
  mutate(
    `Number of Missing Endpoints` = sum(is.na(c_across(unlist(sapply(met_cols,names),use.names = F))))
  ) %>% 
  ungroup() %>% 
  mutate(
    Year_n = (max(Year)-min(Year))-(max(Year)-Year)+1,
    year_class = case_when(
      Year %in% 1995:2000 ~ "1995-2000",
      Year %in% 2001:2006 ~ "2001-2006",
      Year %in% 2007:2012 ~ "2007-2012",
      Year %in% 2013:2018 ~ "2013-2018",
      Year %in% 2019:2024 ~ "2019-2024",
    )
  ) %>% 
  mutate(
    year_class = factor(year_class,levels=rev(c("1995-2000","2001-2006","2007-2012","2013-2018","2019-2024")))
  )


# Filter Data -------------------------------------------------------------
dat <- dat %>% 
  filter(`Years of Data` > 3 & length(unique(year_class)) > 1) %>% 
  group_by(SiteCode,TRCA) %>% 
  filter(
    !TRCA | `Years of Data` > 3,
    !TRCA | any(Year > 2018),
    !TRCA | length(unique(year_class)) > 1,
    !TRCA | grepl("WM$|WMB$",SiteCode)
  ) %>% 
  ungroup() 

# Duplicate TRCA data in full analysis subset -----------------------------
dat <- bind_rows(
  dat %>% mutate(TRCA=F),
  dat %>% filter(TRCA) %>% mutate(TRCA=T)
)

if (F){
  # Get Ranges -------------------------------------------------------------
  
  dat %>% 
    pivot_longer(all_of(unlist(sapply(met_cols,names),use.names = F)),
                 names_to = "Endpoint") %>% 
    ggplot(aes(x=value)) +
    geom_histogram() +
    facet_wrap(~Endpoint,scales="free")
  
  dat %>% 
    pivot_longer(all_of(unlist(sapply(met_cols,names),use.names = F)),
                 names_to = "Endpoint") %>% 
    group_by(Endpoint) %>% 
    summarise(
      min = min(value,na.rm = T),
      max = max(value,na.rm = T)
    ) %>% 
    print(n=Inf)
  
}

# Convert Percent to Proportion -------------------------------------------

dat <- dat %>% 
  mutate(
    across(
      contains(c("Cover","Percent")),
      ~ .x/100
    )
  )

# Replace problematic values ------------------------------------------------


dat <- dat %>% 
  mutate(
    `Habitat Stability Score` = `Habitat Stability Score` - 0.5,
    `Habitat Homogeneity` = 1-`Habitat Homogeneity`,
    across(
      c(D15,D50,D85,D50Max),
      ~ case_when(
        .x < 0 ~ NA_real_,
        .x <= 2 ~ 0,
        .x >= 1000 ~ NA_real_,
        T ~ .x
      )
    ),
    D15 = case_when(
      D15 > 50 ~ NA_real_,
      T ~ D15
    ),
    across(
      c(`Sediment Transport Potential`),
      ~ case_when(
        .x < 0 ~ NA_real_,
        T ~ .x
      )
    ),
    across(
      c(`Sediment Sorting Index`),
      ~ case_when(
        .x < 0 ~ NA_real_,
        T ~ .x
      )
    ),
    across(
      c(
        contains(c("Cover","Percent")),
        `Habitat Stability Score`,
        `Habitat Homogeneity`,
        `Bank Stability - Mean`
      ),
      ~ case_when(
        .x >= 1 ~ quantile(.x[.x<1 & .x>0],0.99,na.rm=T),
        .x < 0 ~ NA_real_,
        T ~ .x
      )
    ),
    across(
      c(
        contains(c("Cover","Percent"))
      ),
      ~ case_when(
        .x < 0.02 ~ 0,
        T ~ .x
      )
    )
  )


# Convert to sf, and attribute with watershed info

dat <- dat %>% 
  sf::st_as_sf(coords=c("Longitude","Latitude"),crs="WGS84",remove = F)


# Geospatial analysis -----------------------------------------------------

TRCA_subwats <- sf::read_sf("data/Subwatersheds_TRCA.geojson")

sf::sf_use_s2(F) # there is an invalid geometry

unzip("data/AEC_V3_Core_Rev0_GPKG_Package02_LakeOntario.zip",exdir=tempdir())
unzip("data/AEC_V3_Core_Rev0_GPKG_Package04_LakeHuronSouth.zip",exdir=tempdir())
unzip("data/AEC_V3_Core_Rev0_GPKG_Package06_LakeSuperior.zip",exdir=tempdir())
unzip("data/AEC_V3_Core_Rev0_GPKG_Package03_OttawaStLawrenceRivers.zip",exdir=tempdir())
strm_lns <- list(
  sf::read_sf(
    file.path(tempdir(),"AEC_V3_Core_Rev0_GPKG_Package02_LakeOntario/w03_Lake_Ontario_West/w03_AEC_V3_Core_Rev0.gpkg"),
    "w03_Reach_SimpleGeometry_Core")
  ,
  sf::read_sf(
    file.path(tempdir(),"AEC_V3_Core_Rev0_GPKG_Package02_LakeOntario/w06_Lake_Ontario_East/w06_AEC_V3_Core_Rev0.gpkg"),
    "w06_Reach_SimpleGeometry_Core")
  ,
  sf::read_sf(
    file.path(tempdir(),"AEC_V3_Core_Rev0_GPKG_Package02_LakeOntario/w04_Lake_Ontario_Central/w04_AEC_V3_Core_Rev0.gpkg"),
    "w04_Reach_SimpleGeometry_Core")
  ,
  sf::read_sf(
    file.path(tempdir(),"AEC_V3_Core_Rev0_GPKG_Package06_LakeSuperior/w23_Lake_Superior_North_West/w23_AEC_V3_Core_Rev0.gpkg"),
    "w23_Reach_SimpleGeometry_Core")
  ,
  sf::read_sf(
    file.path(tempdir(),"AEC_V3_Core_Rev0_GPKG_Package03_OttawaStLawrenceRivers/w08_Ottawa_River_Lower/w08_AEC_V3_Core_Rev0.gpkg"),
    "w08_Reach_SimpleGeometry_Core")
  ,
  sf::read_sf(
    file.path(tempdir(),"AEC_V3_Core_Rev0_GPKG_Package03_OttawaStLawrenceRivers/w09_Ottawa_River_Central_South/w09_AEC_V3_Core_Rev0.gpkg"),
    "w09_Reach_SimpleGeometry_Core")
  ,
  sf::read_sf(
    file.path(tempdir(),"AEC_V3_Core_Rev0_GPKG_Package04_LakeHuronSouth/w14_Georgian_Bay_South_Simcoe/w14_AEC_V3_Core_Rev0.gpkg"),
    "w14_Reach_SimpleGeometry_Core")
  ,
  sf::read_sf(
    file.path(tempdir(),"AEC_V3_Core_Rev0_GPKG_Package04_LakeHuronSouth/w15_Georgian_Bay_Central/w15_AEC_V3_Core_Rev0.gpkg"),
    "w15_Reach_SimpleGeometry_Core")
)

strm_lns <- bind_rows(strm_lns)

strm_lns <- strm_lns %>% 
  filter(Network_Line_Type != "Shoreline Virtual Connector",
         Network_Line_Type != "Virtual Connector") %>% 
  sf::st_transform(sf::st_crs(dat))

dat <- sf::st_join(dat,
                   TRCA_subwats %>% select(WSHED,Subshed)
)

dat <- sf::st_join(dat,
                   strm_lns %>% select(WorkUnitID,ProvReachID,`Strahler Order`=Strahler_Order,`Upstream Catchment Area`=Upstream_Catchment_Area,`Baseflow Index`=BFI_RCA,`Reach Slope`=Slope_ReachChannel_Percent,Turbidity_percUpstreamChannel_T,Turbidity,Flow_Velocity),
                   join =nngeo::st_nn,k=1,maxdist=Inf,parallel=8
)

dat <- dat %>% 
  left_join(
    read_csv("data/AEC_supplementary_MWP_human_disturbance_gradients.csv") %>% 
      select(ProvReachID = ReachID, `Upstream catchment developed (%)`,
             `Upstream catchment agricultural (%)`,
             `Upstream catchment stress index`,
             `Upstream catchment disturbance status`)
  )

dat <- dat %>% 
  mutate(
    WSHED = case_when(
      is.na(WSHED) ~ as.character(WorkUnitID),
      T ~ WSHED
    ),
    Subshed = case_when(
      is.na(Subshed) ~ OrgNameList,
      T ~ Subshed
    )
  )

if (F) mapview::mapview(list(strm_lns,dat))


# Summarize TRCA Streamlines ----------------------------------------------


TRCA_strm_lns <- strm_lns %>% 
  filter(WorkUnitID %in% unique(dat$WorkUnitID)) %>% 
  sf::st_join(TRCA_subwats %>% select(WSHED,Subshed)) %>% 
  left_join(
    read_csv("data/AEC_supplementary_MWP_human_disturbance_gradients.csv") %>% 
      select(ProvReachID = ReachID, `Upstream catchment developed (%)`,
             `Upstream catchment agricultural (%)`,
             `Upstream catchment stress index`,
             `Upstream catchment disturbance status`)
  ) %>% 
  rename(
    `Strahler Order` = Strahler_Order,
    `Upstream Catchment Area`=Upstream_Catchment_Area,
    `Baseflow Index`=BFI_RCA,
    `Reach Slope`=Slope_ReachChannel_Percent
  )

TRCA_strm_lns <- bind_rows(
  TRCA_strm_lns %>% 
    filter(!is.na(Subshed)) %>% 
    mutate(TRCA=T),
  TRCA_strm_lns %>% 
    mutate(TRCA=F) %>% 
    mutate(WSHED = as.character(WorkUnitID),
           Subshed = as.character(WorkUnitID))
)

TRCA_strm_lns <- TRCA_strm_lns %>% 
  mutate(`Stress Class` = cut(`Upstream catchment stress index`,
                              breaks = 3)) 

TRCA_strm_lns <- TRCA_strm_lns %>% 
  mutate(`Stress Class` = cut(`Upstream catchment stress index`,
                              breaks = 3)) 

levels(TRCA_strm_lns$`Stress Class`) <- c("Low","Medium","High")
TRCA_strm_lns$`Stress Class` <- as.character(TRCA_strm_lns$`Stress Class`)
TRCA_strm_lns$`Stress Class`[is.na(TRCA_strm_lns$`Stress Class`)] <- "North"
TRCA_strm_lns$`Stress Class` <- factor(TRCA_strm_lns$`Stress Class`,levels = c("Low","Medium","High","North"))

TRCA_strm_lns <- map(
  c("Subshed","WSHED","Strahler Order","Stress Class","Strahler Order:Stress Class"),
  ~{TRCA_strm_lns %>% 
      as_tibble() %>% 
      mutate(`Strahler Order:Stress Class` = paste(`Strahler Order`,`Stress Class`,sep = ".")) %>% 
      mutate(group=.x) %>% 
      rename_with(.cols=any_of(.x),~paste0("group_val")) %>% 
      mutate(group_val=as.character(group_val)) %>% 
      group_by(TRCA,group,group_val) %>% 
      summarise(
        across(c(
          `Upstream Catchment Area`,
        ),
        ~mean(log10(.x),na.rm=T)^10
        ),
        across(c(
          `Reach Slope`,
        ),
        ~expm1(mean(log1p(.x),na.rm=T))
        ),
        across(c(
          `Baseflow Index`,
          `Upstream catchment developed (%)`,
          `Upstream catchment agricultural (%)`,
          `Upstream catchment stress index`
        ),
        ~median(.x,na.rm=T)
        )
      )
  }
)

TRCA_strm_lns <- bind_rows(TRCA_strm_lns)

write_csv(TRCA_strm_lns,"data/predictor_summaries.csv")

# Convert locations to factors --------------------------------------------

dat <- dat %>% 
  group_by(WSHED) %>% 
  mutate(Longitude_w = median(Longitude),
         Latitude_w = median(Latitude)) %>% 
  group_by(Subshed) %>% 
  mutate(Longitude_s = median(Longitude),
         Latitude_s = median(Latitude)) %>% 
  mutate(
    WSHED = factor(WSHED,levels=unique(WSHED[order(Longitude_w)])),
    Subshed = factor(Subshed,levels=unique(Subshed[order(Longitude_s)])),
    SiteCode = factor(SiteCode,levels=unique(SiteCode[order(Longitude_s)]))
  ) %>% 
  ungroup()


# Bin  streams by cross-sectional size ------------------------------------

dat <- dat %>% 
  mutate(
    Area = `Width - Mean`*`Depth - Mean`
  ) %>% 
  group_by(SiteCode) %>% 
  mutate(Area = median(Area,na.rm=T)) 

dat <- dat %>% 
  group_by(TRCA) %>% 
  mutate(Area_class = cut(`Usable Area`,
                          breaks = scales::log_breaks(5)(`Usable Area`)
  )) %>% 
  ungroup()


dat <- dat %>% 
  mutate(`Stress Class` = cut(`Upstream catchment stress index`,
                              breaks = 3)) 

levels(dat$`Stress Class`) <- c("Low","Medium","High")
dat$`Stress Class` <- as.character(dat$`Stress Class`)
dat$`Stress Class`[is.na(dat$`Stress Class`)] <- "North"
dat$`Stress Class` <- factor(dat$`Stress Class`,levels = c("Low","Medium","High","North"))

dat <- dat %>% 
  group_by(SiteCode) %>% 
  mutate(`Habitat Stability Class` = median(`Habitat Stability Score`,na.rm=T)) %>% 
  ungroup() %>% 
  mutate(`Habitat Stability Class` = cut(`Habitat Stability Class`,
                                         breaks = 2)) 

levels(dat$`Habitat Stability Class`) <- c("Low","High")
dat$`Habitat Stability Class` <- as.character(dat$`Habitat Stability Class`)
dat$`Habitat Stability Class`[is.na(dat$`Habitat Stability Class`)] <- "Missing"
dat$`Habitat Stability Class` <- factor(dat$`Habitat Stability Class`,levels = c("Low","High","Missing"))

dat$`UCA Class` <- cut(dat$`Upstream Catchment Area`,c(0,20,80,Inf))
levels(dat$`UCA Class`) <- c("Small","Intermediate","Large")

# Pivot data long ---------------------------------------------------------

dat_long <- dat %>% 
  pivot_longer(all_of(unlist(sapply(met_cols,names),use.names = F)),
               names_to = "Endpoint") 

# Separate locations ------------------------------------------------------

locs <- dat %>% 
  select(all_of(c(ID_cols[!ID_cols %in% c("ID","SampleDate")],"Sample Years","Years of Data","Sample Months","Sample Year-Months","Longitude","Latitude"))) %>% 
  distinct() %>% 
  sf::st_as_sf(coords=c("Longitude","Latitude"),crs="WGS84")

if (F) mapview::mapview(locs,zcol="OrgNameList")


# Save Data ---------------------------------------------------------------

saveRDS(
  list(
    wide = dat,
    long = dat_long
  ),
  "data/1_formatted_data.RDS"
)

# Calculate Fish Endpoints ------------------------------------------------


lookup_tbl<-read.csv(file.path("data","fish_lookup.csv")) %>% 
  mutate(name.to.use=trimws(name.to.use))
colnames(lookup_tbl)[colnames(lookup_tbl)=="OMNR.Code"]<-"SpeciesCode"

SATI_df<-read.csv(file.path("data","lkpSpeciesHabitatRequirement.csv")) %>% 
  select(SpeciesCode=FWISFishCode,OPT_TEMP) %>% 
  distinct() %>% 
  filter(!is.na(OPT_TEMP))

lookup_tbl<-left_join(lookup_tbl,SATI_df,by="SpeciesCode")

raw_tbl<-read.csv(file.path("data","tblFishSummaryOfTotalCatches.csv")) %>% 
  filter(SiteCode %in% locs$SiteCode,
         year(as.Date(SampleDate))>1995)

raw_tbl<-raw_tbl %>% 
  filter( 
    OSAPSE==1, #the sample event was associated with at least one OSAP project and the sample event itself used site boundaries that were defined as per OSAP
    !is.na(TotalWeightPer100m2) # This makes sure only valid single taxa are included (i.e., Cyprid bulk samples will be excluded)
  )

# Setup Sample Event table ------------------------------------------------

SampleEventID<-raw_tbl[!is.na(raw_tbl$UsableArea) & raw_tbl$UsableArea>0,
                       c("SampleEventID",
                         "StreamName",
                         "StreamCode",
                         "SiteCode",
                         "Latitude",
                         "Longitude",
                         "SampleEventType",
                         "SampleDate")]

SampleEventID<-distinct(SampleEventID)

SampleEventID<-SampleEventID[SampleEventID$SampleEventType=="(Electro) Fish Sampling",]
SampleEventID<-SampleEventID[!is.na(SampleEventID$Latitude) & !is.na(SampleEventID$Longitude),]

# Setup Fish table --------------------------------------------------------

sub_tbl<-raw_tbl[raw_tbl$SampleEventID %in% SampleEventID$SampleEventID,
                 c("SampleEventID","SampleDate","SpeciesCode","CommonName","ScientificName","NumberOfFishPer100m2","TotalWeightPer100m2")] 

# sub_tbl$EstimatedBiomass<-sub_tbl$BulkWeight/sub_tbl$UsableArea*100
# sub_tbl$NumberOfFish<-sub_tbl$NumberOfFish/sub_tbl$UsableArea*100
sub_tbl$EstimatedBiomass<-sub_tbl$TotalWeightPer100m2
sub_tbl$NumberOfFish<-sub_tbl$NumberOfFishPer100m2

sub_tbl<-sub_tbl[sub_tbl$NumberOfFish>0,] %>% 
  group_by(SampleEventID, SampleDate, SpeciesCode,CommonName,ScientificName) %>% 
  summarise(
    NumberOfFish=sum(NumberOfFish,na.rm=T),
    EstimatedBiomass=sum(EstimatedBiomass,na.rm=T)
  ) %>% 
  ungroup() %>% 
  na.omit()

atr_tbl<-left_join(sub_tbl,lookup_tbl,by="SpeciesCode") 



# Taxa Tables -------------------------------------------------------------

taxa_tbl1 <- atr_tbl %>% 
  group_by(SampleEventID,SampleDate) %>% 
  mutate(Comm_Biomass=sum(EstimatedBiomass,na.rm=T),
         Comm_Abundance=sum(NumberOfFish,na.rm=T),
         Perc_Biomass=EstimatedBiomass/Comm_Biomass,
         Perc_Abundance=NumberOfFish/Comm_Abundance
  ) %>% 
  ungroup() %>% 
  select(SampleEventID,SampleDate,SpeciesCode,Comm_Biomass,Comm_Abundance,Perc_Biomass,Perc_Abundance) %>% 
  mutate(SpeciesCode=factor(SpeciesCode)) %>% 
  group_by(SampleEventID,SampleDate) %>% 
  complete(SpeciesCode,fill = list(Comm_Biomass=0,Comm_Abundance=0,Perc_Biomass=0,Perc_Abundance=0)) %>% 
  ungroup() %>% 
  mutate(Comm_Biomass=ifelse(Perc_Abundance>0 & Comm_Biomass==0,NA_real_,Comm_Biomass)) %>% 
  mutate(Perc_Biomass=ifelse(Perc_Abundance>0 & Perc_Biomass==0,NA_real_,Perc_Biomass))


# No Catch Table ----------------------------------------------------------

no_catch_tbl<-SampleEventID %>% 
  as_tibble() %>% 
  select(SampleEventID,SampleDate) %>% 
  left_join(taxa_tbl1,by = c("SampleEventID", "SampleDate")) %>% 
  filter(is.na(Comm_Biomass)) %>% 
  mutate_at(vars(Comm_Biomass:Perc_Abundance),~0) %>% 
  mutate(SpeciesCode=atr_tbl$SpeciesCode[[1]]) %>% 
  mutate(SpeciesCode=factor(SpeciesCode,levels=unique(taxa_tbl1$SpeciesCode))) %>% 
  group_by(SampleEventID,SampleDate) %>% 
  #complete(SpeciesCode,fill = list(Comm_Biomass=0,Comm_Abundance=0,Perc_Biomass=0,Perc_Abundance=0)) %>% 
  ungroup() 

lookup_tbl2<-read_csv(file.path("data","fish_lookup.csv")) %>% 
  mutate(`name to use`=trimws(`name to use`))
colnames(lookup_tbl2)[colnames(lookup_tbl2)=="OMNR Code"]<-"SpeciesCode"

taxa_tbl<-taxa_tbl1 %>% 
  bind_rows(no_catch_tbl) %>% 
  mutate(SpeciesCode=as.character(SpeciesCode)) %>% 
  mutate(SpeciesCode=as.integer(SpeciesCode)) %>% 
  left_join(lookup_tbl2,by="SpeciesCode") %>% 
  group_by(SampleEventID,SampleDate,`name to use`) %>% 
  reframe(across(where(is.numeric),~sum(.,na.rm=T)),
          across(where(is.character),~head(.[!is.na(.)],1))
  ) %>% 
  ungroup() %>% 
  select(any_of(colnames(lookup_tbl2)),everything())

# Calculate endpoints -----------------------------------------------------

atr_tbl<-left_join(sub_tbl,lookup_tbl,by="SpeciesCode") 

ep_tbl<-calc_fish_ep(atr_tbl)

# No Catch Table ----------------------------------------------------------

no_catch_tbl<-SampleEventID %>% 
  select(SampleEventID,SampleDate) %>% 
  left_join(ep_tbl,by = c("SampleEventID", "SampleDate")) %>% 
  filter(is.na(Comm_Richness)) %>% 
  mutate_at(vars(Comm_Biomass:Rep_Lith_Rich),~0)

ep_tbl[sapply(ep_tbl,is.nan)]<-NA

ep_tbl<-ep_tbl %>% 
  bind_rows(no_catch_tbl)


# Endpoint Index ----------------------------------------------------------
ep_tbl_lng<-ep_tbl %>% 
  gather("Endpoint","Value",-SampleEventID,-SampleDate) %>% 
  mutate(Value=ifelse(is.nan(Value),NA,Value)) %>% 
  mutate(Endpoint_group=case_when(
    grepl("^Therm_",Endpoint) ~ "Thermal Regime",
    grepl("^Troph_",Endpoint) ~ "Trophic Guild",
    grepl("^Tol_",Endpoint) ~ "Tolerance Group",
    grepl("^Spec_",Endpoint) ~ "Individual Species",
    grepl("^Env_",Endpoint) ~ "Environment Preference",
    grepl("^Rep_",Endpoint) ~ "Reproductive Guild",
    grepl("^Comm_",Endpoint) ~ "General Community"
  )) %>%
  mutate(Endpoint_subgroup=str_split(Endpoint,"_",simplify = T)[,2]) %>% 
  mutate(Endpoint_subgroup=case_when(
    Endpoint_subgroup == "SATI" ~ "SATI",
    Endpoint_subgroup == "cold" ~ "Cold water",
    Endpoint_subgroup == "cool" ~ "Cool water",
    Endpoint_subgroup == "coolcold" ~ "Cold-Cool water",
    Endpoint_subgroup == "warm" ~ "Warm water",
    Endpoint_subgroup == "invert" ~ "Invertivore",
    Endpoint_subgroup == "detrit" ~ "Detritivore",
    Endpoint_subgroup == "carn" ~ "Carniivore",
    Endpoint_subgroup == "herb" ~ "Herbiivore",
    Endpoint_subgroup == "plank" ~ "Planktivore",
    Endpoint_subgroup == "tol" ~ "Tolerant",
    Endpoint_subgroup == "intol" ~ "Intolerant",
    Endpoint_subgroup == "RBD" ~ "Rainbow Darter",
    Endpoint_subgroup == "BKT" ~ "Brook Trout",
    Endpoint_subgroup == "GDF" ~ "Goldfish",
    Endpoint_subgroup == "RDG" ~ "Round Goby",
    Endpoint_subgroup == "CKC" ~ "Creek Chub",
    Endpoint_subgroup == "CMC" ~ "Common Carp",
    Endpoint_subgroup == "SKP" ~ "Mottled and Slimy Sculpin",
    Endpoint_subgroup == "LMP" ~ "Lampreys",
    Endpoint_subgroup == "DRT" ~ "Sensitive Darters",
    Endpoint_subgroup == "CEB" ~ "Centrarchids and Basses",
    Endpoint_subgroup == "CTF" ~ "Catfishes",
    Endpoint_subgroup == "SAL" ~ "Salmonids",
    Endpoint_subgroup == "benthic" ~ "Benthic",
    Endpoint_subgroup == "pelagic" ~ "Pelagic",
    Endpoint_subgroup == "benthopelagic" ~ "Benthopelagic",
    Endpoint_subgroup == "guard" ~ "Guarders",
    Endpoint_subgroup == "BrHid" ~ "Brood Hiders/Nest Spawners",
    Endpoint_subgroup == "OSS" ~ "Open Substratum/Substratum Choosers",
    Endpoint_subgroup == "Lith" ~ "Lithophilic Spawners",
    Endpoint=="Comm_Biomass" ~ "Total Biomass",
    Endpoint=="Comm_Abundance" ~ "Total Abundance",
    Endpoint=="Comm_Richness" ~ "Total Richness"
  )) %>% 
  mutate(Endpoint_type=case_when(
    grepl("_Bioperc$",Endpoint) ~ "Biomass (%)",
    grepl("_Abuperc$",Endpoint) ~ "Abundance (%)",
    grepl("_Bio$",Endpoint) ~ "Total Biomass",
    grepl("_Abu$",Endpoint) ~ "Total Abundance",
    grepl("_Rich$",Endpoint) ~ "Richness",
    Endpoint=="Comm_Biomass" ~ "Total Biomass",
    Endpoint=="Comm_Abundance" ~ "Total Abundance",
    Endpoint=="Comm_Richness" ~ "Richness"
  )) %>% 
  filter(!is.na(Value))


# Write results to excel --------------------------------------------------

ep_tbl_out<-ep_tbl %>% 
  left_join(SampleEventID) %>% 
  select(StreamName, StreamCode, SiteCode, Latitude, Longitude, SampleEventType, SampleDate,SampleEventID,everything())

if (F){
  ep_tbl_out %>% 
    mutate(Year = year(as.Date(SampleDate)),
           month = month(as.Date(SampleDate))) %>% 
    group_by(SiteCode,Year) %>% 
    summarise(
      n_samples = n(),
      months = paste0(month,collapse = ", ")
    ) %>% 
    arrange(desc(n_samples)) %>% 
    View()
}

ep_tbl_lng_out<-ep_tbl_lng %>% 
  left_join(SampleEventID) %>% 
  select(StreamName, StreamCode, SiteCode, Latitude, Longitude, SampleEventType, SampleDate,SampleEventID,everything())

param_tbl_out<-ep_tbl_lng %>% 
  select(Endpoint_group,Endpoint_subgroup,Endpoint_type,Endpoint) %>% 
  distinct()

taxa_tbl_out<-taxa_tbl %>% 
  left_join(SampleEventID) %>% 
  select(StreamName, StreamCode, SiteCode, Latitude, Longitude, SampleEventType, SampleDate,SampleEventID,everything())


wb = createWorkbook()

addWorksheet(wb, "Endpoints")
addWorksheet(wb, "Results - Wide")
addWorksheet(wb, "Results - Long")
addWorksheet(wb, "Taxa Table")

writeData(wb,param_tbl_out, sheet="Endpoints", rowNames =FALSE)
writeData(wb,ep_tbl_out, sheet="Results - Wide", rowNames =FALSE)
writeData(wb,ep_tbl_lng_out, sheet="Results - Long", rowNames =FALSE)
writeData(wb,taxa_tbl_out, sheet="Taxa Table", rowNames =FALSE)

saveWorkbook(wb, file.path("data","Fish_Endpoints.xlsx"),overwrite = T)

final_fish_out <- ep_tbl_lng_out %>% 
  mutate(Year = year(as.Date(SampleDate))) %>% 
  select(SiteCode,Year,Endpoint,Endpoint_group,Endpoint_subgroup,Endpoint_type,Value) %>% 
  group_by(SiteCode,Year,Endpoint,Endpoint_group,Endpoint_subgroup,Endpoint_type) %>% 
  summarise(Value = median(Value,na.rm = T), .groups="drop") 

write_csv(final_fish_out,"data/fish_endpoints.csv")

4 Fit Models and Format Results

source("00_source.R")
library(glmmTMB)
library(DHARMa)
library(broom.mixed)
library(parameters)

dat_in <- readRDS("data/1_formatted_data.RDS")$long 
predictor_summaries<-read_csv("data/predictor_summaries.csv")

if (F) {
  ggplot(dat_in,aes(x=value)) +
    geom_histogram() +
    facet_wrap(~Endpoint,scales="free")
  
  ggplot(dat_in,aes(x=value)) +
    geom_histogram() +
    scale_x_continuous(trans="log1p")+
    facet_wrap(~Endpoint,scales="free")
}

# Setup Model param -------------------------------------------------------

dat_in <- dat_in %>% 
  filter(`Stress Class` != "North") %>% # Exclude sites without land cover data
  #filter(Endpoint != c("Habitat Stability Score")) %>%  # Has a problematic distribution
  filter(Endpoint != c("Sediment Transport Potential")) %>%  # Has a problematic distribution
  group_by(TRCA) %>% 
  mutate(
    Year_n2 = Year,
  ) %>% 
  mutate(
    Year_n = (max(Year)-min(Year))-(max(Year)-Year)+1,
  ) %>% 
  ungroup() %>% 
  mutate(
    `Strahler Order` = as.factor(`Strahler Order`)
  )

ep_tbl <- enframe(met_cols) %>% 
  unnest(value) %>% 
  mutate(Endpoint = names(value)) %>% 
  unnest(value)

mod_data_fin <- dat_in %>% 
  as_tibble() %>% 
  group_by(across(any_of(c("TRCA","Family","Endpoint")))) %>% 
  nest() %>% 
  ungroup() %>% 
  left_join(enframe(met_cols) %>% 
              unnest(value) %>% 
              mutate(Endpoint = names(value)) %>% 
              unnest(value) %>% 
              select(-value,Endpoint_group=name,Endpoint),
            by = "Endpoint" ) %>% 
  mutate(Endpoint = factor(Endpoint,levels = ep_tbl$Endpoint),
         Endpoint_group = factor(Endpoint_group,levels = unique(ep_tbl$name))) %>% 
  mutate(model_params = case_when(
    Endpoint %in% c("D15",
                    "D50",
                    "D85",
                    "D50Max"
    ) ~ list(
      list(
        family = lognormal(link = "log"),
        ziformula = ~ 1,
        dispformula = ~ 1
      )
    ),
    Endpoint %in% c("Sediment Sorting Index"
    ) ~ list(
      list(
        family = ziGamma(link = "log"),
        ziformula = ~ 1,
        dispformula = ~ 1
      )
    ),
    Endpoint %in% c("Width - Mean",
                    "Depth - Mean") ~ list(
                      list(
                        family = lognormal(link = "log"),
                        ziformula = ~ 0,
                        dispformula = ~ 1
                      )
                    ),
    Endpoint %in% c("Sediment Transport Potential") ~ list(
      list(
        family = Gamma(link = "sqrt"),
        ziformula = ~ 0,
        dispformula = ~ 1
      )
    ),
    Endpoint %in% c("Usable Area","Width:Depth Ratio") ~ list(
      list(
        family = Gamma(link = "log"),
        ziformula = ~ 0,
        dispformula = ~ 1
      )
    ),
    grepl(" Cover| Percent",Endpoint) | Endpoint %in% c("Habitat Homogeneity",
                                                        "Habitat Stability Score")  ~ list(
                                                          list(
                                                            family = beta_family(link = "logit"),
                                                            ziformula = ~ 1,
                                                            dispformula = ~ 1
                                                          )
                                                        ),
    Endpoint %in% c("Bank Stability - Mean")  ~ list(
      list(
        family = beta_family(link = "probit"),
        ziformula = ~ 0,
        dispformula = ~ 1
      )
    ),
    Endpoint %in% c("") ~ list(
      list(
        family = gaussian(link = "identity"),
        ziformula = ~ 1,
        dispformula = ~ 1
      )
    )
  )) %>% 
  mutate(model_params = map2(TRCA,model_params,
                             ~case_when(.x==T ~ c(.y,
                                                  list(
                                                    control = glmmTMBControl(parallel = parallel::detectCores()),
                                                    formula = as.formula("value ~ log1p(`Reach Slope`) + Year_n*`UCA Class`*`Stress Class` + (Year_n|`UCA Class`:`Stress Class`:SiteCode) + (1|Year_n)")
                                                  )),
                                        .x==F ~ c(.y,
                                                  list(
                                                    control = glmmTMBControl(parallel = parallel::detectCores()),
                                                    formula = as.formula("value ~ log1p(`Reach Slope`) + Year_n*`UCA Class`*`Stress Class` + (Year_n|`UCA Class`:`Stress Class`:SiteCode) + (1|Year_n)")
                                                  ))
                             )
                             
  )) %>% 
  mutate(model_params = map2(TRCA,model_params,
                             ~{
                               if (.y$dispformula[[2]]!=0 & .x==T){
                                 .y$dispformula <- as.formula(" ~ 1")
                               }
                               if (.y$dispformula[[2]]!=0 & .x==F){
                                 .y$dispformula <- as.formula(" ~ 1")
                               }
                               if (.y$ziformula[[2]]!=0 & .x==T){
                                 .y$ziformula <- as.formula(" ~ log1p(`Reach Slope`) + Year_n + `UCA Class` + `Stress Class` + (Year_n|`UCA Class`:`Stress Class`:SiteCode) + (1|Year_n)")
                               }
                               if (.y$ziformula[[2]]!=0 & .x==F){
                                 .y$ziformula <- as.formula(" ~ log1p(`Reach Slope`) + Year_n + `UCA Class` + `Stress Class` + (Year_n|`UCA Class`:`Stress Class`:SiteCode) + (1|Year_n)")
                               }
                               return(.y)
                             }
  )) %>% 
  mutate(model_params = map2(model_params,data, 
                             ~c(.x,
                                list(
                                  data= .y %>% 
                                    filter(!is.na(value)) %>% 
                                    select(
                                      any_of(
                                        c("WSHED","Subshed","SiteCode","Area_class",
                                          "Sample_Season","Year_n","Year_n2","value",
                                          "year_class","Year",
                                          "Upstream Catchment Area","Baseflow Index","Reach Slope","Strahler Order",
                                          "Upstream catchment developed (%)","Upstream catchment agricultural (%)",
                                          "Upstream catchment stress index","Stress Class","Habitat Stability Class","UCA Class")
                                      )
                                    ) 
                                )
                             ))) %>% 
  mutate(Endpoint = factor(Endpoint, levels = unlist(sapply(met_cols,names),use.names=F)))


mod_data_fin <- mod_data_fin %>% 
  mutate(p0 = map_dbl(data,~sum(.x$value==0,na.rm=T)/length(.x$value==0)))

mod_data_fin <- mod_data_fin %>% 
  filter(p0 < 0.75)

# Run Models --------------------------------------------------------------

mod_data_fin <- mod_data_fin %>% 
  mutate(
    model = map(model_params, ~do.call(glmmTMB::glmmTMB,.x))
  ) 

mod_data_fin <- mod_data_fin %>% 
  mutate(
    pred_error = map_lgl(model,
                         ~any(
                           is.na(
                             predict(.x,
                                     type="link",
                                     re.form = NULL,
                                     se.fit = TRUE,
                                     cov.fit = F,
                                     do.bias.correct = F,
                                     allow.new.levels = T)$se.fit
                           ) |
                             any(is.na(as.data.frame(glmmTMB::ranef(.x, condVar = T))$condsd))
                         ) 
    )
  )

print(mod_data_fin,n=Inf)

# Exclude outliers --------------------------------------------------------
mod_data_fin <- mod_data_fin %>% 
  mutate(
    model_params = case_when(
      pred_error == T ~ map2(model_params,model, 
                             ~{
                               .x$data <- .x$data[!DHARMa::outliers(DHARMa::simulateResiduals(.y, use.u = T,refit = F,plot = F,n=9999),return = "logical",lowerQuantile = 0.025, upperQuantile = 0.975),]
                               return(.x)
                             }
      ),
      T ~ model_params
    )
  ) 

# Refit Models --------------------------------------------------------------

mod_data_fin <- mod_data_fin %>%
  mutate(
    model = case_when(
      pred_error == T ~ map(model_params, ~do.call(glmmTMB::glmmTMB,.x)), #refit models with default links
      T ~ model
    )
  )

mod_data_fin <- mod_data_fin %>% 
  mutate(
    pred_error = map_lgl(model,
                         ~any(
                           is.na(
                             predict(.x,
                                     type="link",
                                     re.form = NULL,
                                     se.fit = TRUE,
                                     cov.fit = F,
                                     do.bias.correct = F,
                                     allow.new.levels = T)$se.fit
                           ) |
                             any(is.na(as.data.frame(glmmTMB::ranef(.x, condVar = T))$condsd))
                         ) 
    )
  )



print(mod_data_fin,n=Inf)

# adjust link functions for bad fits
mod_data_fin <- mod_data_fin %>% 
  mutate(model = pmap(
    list(.x=model,
         .y=pred_error,
         .z=model_params),
    function(.x,.y,.z){
      if (!.y) return(.x)
      #if (.x$call$family$family != "beta") return(.x)
      if (.x$call$family$family == "beta") {
        .z$family <- beta_family(link = "probit")
        .x <- do.call(glmmTMB::glmmTMB,.z)
      }
      if (.x$call$family$family == "Gamma") {
        .z$family <- ziGamma(link = "sqrt")
        .x <- do.call(glmmTMB::glmmTMB,.z)
      }
      if (.x$call$family$family == "lognormal") {
        .z$family <- lognormal(link = "sqrt")
        .x <- do.call(glmmTMB::glmmTMB,.z)
      }
      return(.x)
    })
  )

mod_data_fin <- mod_data_fin %>% 
  mutate(
    pred_error = map2_lgl(model,
                          pred_error,
                          ~{
                            if (!.y) return(.y)
                            any(
                              is.na(
                                predict(.x,
                                        type="link",
                                        re.form = NULL,
                                        se.fit = TRUE,
                                        cov.fit = F,
                                        do.bias.correct = F,
                                        allow.new.levels = T)$se.fit
                              )
                            ) |
                              any(is.na(as.data.frame(glmmTMB::ranef(.x, condVar = T))$condsd))
                          })
  )

print(mod_data_fin,n=Inf)



mod_data_fin <- mod_data_fin %>% 
  filter(!pred_error)

# Predict observations ----------------------------------------------------

mod_data_fin <- mod_data_fin %>% 
  mutate(Pred1 = map2(model,model_params, function (i,nm) {
    
    year_vals <- nm$data %>% 
      select(Year_n,Year_n2,year_class) %>% 
      distinct()  
    
    pred_frame <- nm$data %>% 
      select(any_of(c("WSHED","Subshed","Area_class","SiteCode",
                      "Upstream Catchment Area","Baseflow Index","Reach Slope","Strahler Order",
                      "Upstream catchment developed (%)","Upstream catchment agricultural (%)",
                      "Upstream catchment stress index","Stress Class","Habitat Stability Class","UCA Class"))) %>% 
      distinct() %>% 
      mutate(Sample_Season="Summer",#list(factor(levels(nm$data$Sample_Season),levels = levels(nm$data$Sample_Season))),
             Year_n = list(seq(from=min(year_vals$Year_n),to=max(year_vals$Year_n)))
      ) %>%
      unnest(Year_n) %>%
      unnest(Sample_Season) %>% 
      left_join(
        year_vals
      )
    
    pred_frame<-pred_frame[complete.cases(pred_frame),]
    
    a1 <- pred_frame %>% 
      bind_cols(
        data.frame(
          pred=predict(i,
                       newdata = pred_frame,
                       type="link",
                       re.form = NULL,
                       se.fit = TRUE,
                       cov.fit = F,
                       do.bias.correct = F,
                       allow.new.levels = T)
        )
      ) %>%
      mutate(pred.fit = predict(i,newdata = pred_frame,type="response")) %>% 
      mutate(pred.fit = i$modelInfo$family$linkfun(pred.fit)) %>% 
      mutate(
        pred_low = pred.fit - (qnorm(1 - 0.1 / 2) * pred.se.fit),
        pred_high = pred.fit + (qnorm(1 - 0.1 / 2) * pred.se.fit),
        pred_low_trend = pred.fit - (qnorm(1 - 0.5 / 2) * pred.se.fit),
        pred_high_trend = pred.fit + (qnorm(1 - 0.5 / 2) * pred.se.fit),
      ) %>% 
      left_join(
        nm$data %>% select(any_of(c("WSHED","Subshed","SiteCode","Area_class",
                                    "Sample_Season","Year_n","Year_n2",
                                    "value","year_class",
                                    "Upstream Catchment Area","Baseflow Index","Reach Slope","Strahler Order",
                                    "Upstream catchment developed (%)","Upstream catchment agricultural (%)",
                                    "Upstream catchment stress index","Stress Class","Habitat Stability Class","UCA Class")))
      ) %>% 
      mutate(
        across(
          c(pred.fit,pred_low,pred_high,pred_low_trend,pred_high_trend),
          ~i$modelInfo$family$linkinv(.x)
        )
      ) 
    
    return(a1)
  }))

mod_data_fin <- mod_data_fin %>% 
  mutate(Pred1_error = map_lgl(Pred1,~any(is.na(.x$pred.se.fit))))

print(mod_data_fin,n=Inf)

# Extract random effects terms --------------------------------------------

mod_data_fin <- mod_data_fin %>% 
  mutate(rand_ef = map(model,~get_terms2(.x))) 


mod_data_fin <- mod_data_fin %>% 
  mutate(rand_ef_error = map_lgl(rand_ef,~any(is.na(.x$sd))))

print(mod_data_fin,n=Inf)

# Get simulated residuals -------------------------------------------------

mod_data_fin <- mod_data_fin %>%
  mutate(
    sim_resid_refitF = map(model, ~ DHARMa::simulateResiduals(.x, use.u = T,refit = F,plot = F,n=9999)),
  )

# Format final data tables ------------------------------------------------

mod_data_fin <- mod_data_fin %>%
  mutate(
    Pred2 = map(Pred1,
                ~ .x %>% 
                  group_by(WSHED, Subshed, Sample_Season, SiteCode, year_class,
                           `Upstream Catchment Area`,`Baseflow Index`,`Reach Slope`,`Strahler Order`,
                           `Upstream catchment developed (%)`,`Upstream catchment agricultural (%)`,
                           `Upstream catchment stress index`,`Stress Class`,`UCA Class`) %>% #,`Habitat Stability Class`
                  summarise(
                    across(
                      c(pred.fit, pred.se.fit, pred_low, pred_high, pred_low_trend, pred_high_trend, value),
                      ~median(.x, na.rm=T)
                    ),
                    .groups = "drop"
                  ) %>% 
                  mutate(`UCA Class`=as.character(`UCA Class`))
    ),
    Pred2 = map2(Pred2,rand_ef,
                 function(x,y){
                   x %>% 
                     left_join(
                       y %>% 
                         filter(term == "Year_n") %>% 
                         filter(group_var== "`UCA Class`:`Stress Class`:SiteCode") %>% #:`Habitat Stability Class`
                         filter(component == "cond") %>% 
                         select(-group_var,-term,-sd,-`1`) %>% 
                         rename(lower_slope_sig = pred_low,
                                upper_slope_sig = pred_high,
                                lower_slope_trend = pred_low_trend,
                                upper_slope_trend = pred_high_trend,
                                mean_slope = value,
                         )
                     )
                 }),
    Pred2 = map(Pred2,
                ~.x %>% 
                  mutate(
                    linear_trend = case_when(
                      lower_slope_sig > 0 ~ "Increasing Significant",
                      upper_slope_sig < 0 ~ "Decreasing Significant",
                      lower_slope_trend > 0 ~ "Increasing Trend",
                      upper_slope_trend < 0 ~ "Decreasing Trend",
                      T ~ "No Change"
                    )
                  ) %>% 
                  mutate(
                    linear_trend = factor(linear_trend,
                                          levels = c(
                                            "Increasing Significant",
                                            "Increasing Trend",
                                            "No Change",
                                            "Decreasing Trend",
                                            "Decreasing Significant"
                                          ))
                  )
                
    )
  )


# UCA, Stability and stress class trends ------------------------------------

mod_data_fin <- mod_data_fin %>% 
  mutate(
    spatial_summ = map(data,~summary_fun(.x))
  ) %>% 
  mutate(summary_terms_temporal = map(rand_ef,
                                      ~{
                                        .x %>% 
                                          filter(term == "Year_n") %>% 
                                          filter(group_var!= "`UCA Class`:`Stress Class`:SiteCode") %>% #:`Habitat Stability Class`
                                          filter(component == "fixed") %>% 
                                          rename(lower_slope_sig = pred_low,
                                                 upper_slope_sig = pred_high,
                                                 lower_slope_trend = pred_low_trend,
                                                 upper_slope_trend = pred_high_trend,
                                                 mean_slope = value,
                                          ) %>% 
                                          mutate(
                                            linear_trend = case_when(
                                              lower_slope_sig > 0 ~ "Increasing Significant",
                                              upper_slope_sig < 0 ~ "Decreasing Significant",
                                              lower_slope_trend > 0 ~ "Increasing Trend",
                                              upper_slope_trend < 0 ~ "Decreasing Trend",
                                              T ~ "No Change"
                                            )
                                          ) %>% 
                                          mutate(
                                            linear_trend = factor(linear_trend,
                                                                  levels = c(
                                                                    "Increasing Significant",
                                                                    "Increasing Trend",
                                                                    "No Change",
                                                                    "Decreasing Trend",
                                                                    "Decreasing Significant"
                                                                  ))
                                          ) 
                                      })) %>% 
  mutate(summary_terms_spatial = map(rand_ef,
                                     ~{.x %>% 
                                         filter(term == "(Intercept)") %>% #labelled intercept but is predicted at max year
                                         filter(group_var != "`UCA Class`:`Stress Class`:SiteCode") %>% #:`Habitat Stability Class`
                                         filter(component == "fixed") 
                                     }))


saveRDS(mod_data_fin,"data/results_out.rds")

# Model Fit tables --------------------------------------------------------------

mod_tbl_list <- list()

mod_tbl_list$TRCA_mod_tbl <- modelsummary::modelsummary(
  setNames(
    mod_data_fin$model[mod_data_fin$TRCA],
    mod_data_fin$Endpoint[mod_data_fin$TRCA]
  ),
  vcov = sandwich::vcovHC,
  shape = effect:component + term + group ~ model + statistic,
  stars = TRUE,
  gof_omit = c("BIC|AIC|Std|Cond|Marg|ICC"),
  gof_function = function(model) tibble(`Pseudo R2`=scales::label_number(0.001)(cor(predict(model),model$frame$value,method="spearman")^2)) ,
  output = "gt",
  coef_rename = F,
  notes = "zi = zero-inflation; cond = conditional; fixed = fixed effects parameter; ran_pars = random effects parameters;",
  title = "TRCA Data Only"
) 

sig_list <- which(apply(mod_tbl_list$TRCA_mod_tbl$`_data`,1,function(x)grepl("\\*|\\+",x)),arr.ind = T)
sig_list <- split(sig_list,1:nrow(sig_list))

for (i in sig_list){
  mod_tbl_list$TRCA_mod_tbl <- mod_tbl_list$TRCA_mod_tbl %>% 
    gt::tab_style(
      style = gt::cell_fill(color = "lightblue"),
      locations = gt::cells_body(
        columns = i[[1]],
        rows = i[[2]]
      )
    )
}

if (T){
  mod_tbl_list$ALL_mod_tbl <- modelsummary::modelsummary(
    setNames(
      mod_data_fin$model[!mod_data_fin$TRCA],
      mod_data_fin$Endpoint[!mod_data_fin$TRCA]
    ),
    vcov = sandwich::vcovHC,
    shape = effect:component + term + group ~ model + statistic,
    stars = TRUE,
    gof_omit = c("BIC|AIC|Std|Cond|Marg|ICC"),
    gof_function = function(model) tibble(`Pseudo R2`=scales::label_number(0.001)(cor(predict(model),model$frame$value,method="spearman")^2)),
    output = "gt",
    coef_rename = F,
    notes = "zi = zero-inflation; cond = conditional; fixed = fixed effects parameter; ran_pars = random effects parameters;",
    title = "All Data"
  ) 
  
  sig_list <- which(apply(mod_tbl_list$ALL_mod_tbl$`_data`,1,function(x)grepl("\\*|\\+",x)),arr.ind = T)
  sig_list <- split(sig_list,1:nrow(sig_list))
  
  for (i in sig_list){
    mod_tbl_list$ALL_mod_tbl <- mod_tbl_list$ALL_mod_tbl %>% 
      gt::tab_style(
        style = gt::cell_fill(color = "lightblue"),
        locations = gt::cells_body(
          columns = i[[1]],
          rows = i[[2]]
        )
      )
  }
}


saveRDS(mod_tbl_list,"data/mod_tbl_list.rds")

# Correlation Figures -----------------------------------------------------

cor_list <- list(
  TRCA_spat = mod_data_fin %>% 
    filter(TRCA) %>% 
    select(Endpoint, data) %>% 
    unnest(data) %>% 
    pivot_wider(names_from = Endpoint,values_from = value) %>% 
    select(any_of(unlist(sapply(met_cols,names),use.names=F)),
           any_of(c("Upstream Catchment Area","Reach Slope","Upstream catchment developed (%)","Upstream catchment agricultural (%)","Upstream catchment stress index"))
    ) 
  ,
  All_spat = mod_data_fin %>% 
    filter(!TRCA) %>% 
    select(Endpoint, data) %>% 
    unnest(data) %>% 
    pivot_wider(names_from = Endpoint,values_from = value) %>% 
    select(any_of(unlist(sapply(met_cols,names),use.names=F)),
           any_of(c("Upstream Catchment Area","Reach Slope","Upstream catchment developed (%)","Upstream catchment agricultural (%)","Upstream catchment stress index"))
    )
)

cor_list <- map(cor_list,function(x){
  heatmap_fn(x,
             alpha=0.1,
             sig_est=0.6,
             method="spearman",
             scale_grad=ggplot2::scale_fill_gradient2(low = scales::muted(rgb(red=0/255,green=102/255,blue=255/255)),
                                                      mid = "white",
                                                      high = scales::muted(rgb(red=0/255,green=176/255,blue=80/255)),
                                                      limits=c(-1,1))) 
})



trend_cor <- mod_data_fin %>% 
  filter(TRCA) %>% 
  select(Endpoint, Pred2) %>% 
  mutate(Pred_current = map(Pred2,
                            ~.x %>% 
                              filter(year_class==max(levels(year_class)),
                                     Sample_Season=="Summer"))) %>% 
  select(-Pred2) %>% 
  unnest(Pred_current) %>% 
  select(SiteCode,Endpoint,linear_trend ) %>% 
  mutate(linear_trend = case_when(
    linear_trend == "Increasing Significant" ~ 2,
    linear_trend == "Increasing Trend" ~ 1,
    linear_trend == "No Change" ~ 0,
    linear_trend == "Decreasing Trend" ~ -1,
    linear_trend == "Decreasing Significant" ~ -2,
  )) %>% 
  pivot_wider(values_from = linear_trend ,names_from = Endpoint) %>% 
  select(-SiteCode)

t1 <- cor(trend_cor,method = "spearman",use="pairwise.complete.obs")
keep <- (!apply(t1,2,function(x) all(is.na(x))))
trend_cor <- trend_cor[,keep]

trend_cor<-heatmap_fn(
  trend_cor,
  alpha=0.1,
  sig_est=0.6,
  method="spearman",
  use="pairwise.complete.obs",
  scale_grad=ggplot2::scale_fill_gradient2(low = scales::muted(rgb(red=0/255,green=102/255,blue=255/255)),
                                           mid = "white",
                                           high = scales::muted(rgb(red=0/255,green=176/255,blue=80/255)),
                                           limits=c(-1,1))) 

cor_list$TRCA_temp <- trend_cor

trend_cor <- mod_data_fin %>% 
  filter(!TRCA) %>% 
  select(Endpoint, Pred2) %>% 
  mutate(Pred_current = map(Pred2,
                            ~.x %>% 
                              filter(year_class==max(levels(year_class)),
                                     Sample_Season=="Summer"))) %>% 
  select(-Pred2) %>% 
  unnest(Pred_current) %>% 
  select(SiteCode,Endpoint,linear_trend ) %>% 
  mutate(linear_trend = case_when(
    linear_trend == "Increasing Significant" ~ 2,
    linear_trend == "Increasing Trend" ~ 1,
    linear_trend == "No Change" ~ 0,
    linear_trend == "Decreasing Trend" ~ -1,
    linear_trend == "Decreasing Significant" ~ -2,
  )) %>% 
  pivot_wider(values_from = linear_trend ,names_from = Endpoint) %>% 
  select(-SiteCode)

t1 <- cor(trend_cor,method = "spearman",use="pairwise.complete.obs")
keep <- (!apply(t1,2,function(x) all(is.na(x))))
trend_cor <- trend_cor[,keep]

trend_cor<-heatmap_fn(
  trend_cor,
  alpha=0.1,
  sig_est=0.6,
  method="spearman",
  use="pairwise.complete.obs",
  scale_grad=ggplot2::scale_fill_gradient2(low = scales::muted(rgb(red=0/255,green=102/255,blue=255/255)),
                                           mid = "white",
                                           high = scales::muted(rgb(red=0/255,green=176/255,blue=80/255)),
                                           limits=c(-1,1))) 

cor_list$All_temp <- trend_cor

saveRDS(cor_list,"data/cor_plots.rds")

5 Generate Figures

source("00_source.R")
library(mapview)
library(basemaps)
library(ggspatial)

cols_list <- tail(scales::col_darker(RColorBrewer::brewer.pal(4,"Set1"),5),3)
names(cols_list)<-c("Small","Intermediate","Large")

mod_data_fin <- readRDS("data/results_out.rds")
coords <- readRDS("data/1_formatted_data.RDS")$wide %>% 
  select(TRCA,WSHED,SiteCode,`Stress Class`,`UCA Class`) %>% #,`Habitat Stability Class`
  group_by(TRCA,WSHED,SiteCode,`Stress Class`,`UCA Class`) %>% #,`Habitat Stability Class`
  summarise(geometry = unique(st_centroid(geometry)),.groups="drop")

saveRDS(coords,"data/coords.rds")

coords<-coords %>%
  select(-TRCA,-WSHED,-`Stress Class`,-`UCA Class`) %>% #,-`Habitat Stability Class`
  distinct()

sf::sf_use_s2(F)

TRCA_subwats <- sf::read_sf("data/Subwatersheds_TRCA.geojson") 
TRCA_wats <- TRCA_subwats %>% 
  group_by(WSHED) %>% 
  summarise(geometry = sf::st_union(geometry))

unzip("data/AEC_V3_Core_Rev0_GPKG_Package02_LakeOntario.zip",exdir=tempdir())
strm_lns <- sf::read_sf(
  file.path(tempdir(),"AEC_V3_Core_Rev0_GPKG_Package02_LakeOntario/w03_Lake_Ontario_West/w03_AEC_V3_Core_Rev0.gpkg"),
  "w03_Reach_SimpleGeometry_Core") %>% 
  sf::st_transform(sf::st_crs(TRCA_subwats)) 

strm_lns <- strm_lns[apply(sf::st_intersects(strm_lns,sf::st_buffer(TRCA_wats,-0.0001),sparse=F),1,any),]

strm_lns <- strm_lns %>% 
  filter(Network_Line_Type != "Shoreline Virtual Connector")

mod_data_fin <- mod_data_fin %>% 
  filter(!TRCA)

plot_out <- mod_data_fin %>% 
  mutate(fam = map_chr(model_params,~.x$family$family)) %>% 
  mutate(trans = case_when(
    fam == "beta" ~ list(clogit),
    fam == "lognormal" ~ list(scales::transform_pseudo_log()),
    fam == "Gamma" ~ list(scales::transform_pseudo_log()),
  )) %>% 
  mutate(brks = case_when(
    fam == "beta" ~ list(scales::pretty_breaks),
    fam == "lognormal" ~ list(pseudo_log_breaks),
    fam == "Gamma" ~ list(pseudo_log_breaks),
  )) %>% 
  select(TRCA, Endpoint, Endpoint_group, fam, trans, brks, Pred1, Pred2, summary_terms_temporal, summary_terms_spatial) %>% 
  mutate(Pred_current = map(Pred2,
                            ~left_join(.x,coords,by = join_by(SiteCode)) %>% 
                              filter(year_class==max(levels(year_class)),
                                     Sample_Season=="Summer") %>% 
                              st_as_sf() %>% 
                              filter(is.na(as.numeric(as.character(WSHED)))) %>% 
                              mutate(`UCA Class` = factor(`UCA Class`,levels=c("Small","Intermediate","Large")))
  ))%>% 
  mutate(Pred_obs = map(Pred1,
                        ~.x %>% 
                          distinct() %>% 
                          filter(!is.na(value)) %>% 
                          mutate(Year_n = 2000+Year_n) %>% 
                          filter(is.na(as.numeric(as.character(WSHED)))) %>% 
                          mutate(`UCA Class` = factor(`UCA Class`,levels=c("Small","Intermediate","Large")))
  )) %>% 
  select(-Pred2,-Pred1) 


# Raw Data Plots ----------------------------------------------------------

plot_out$raw_plots <- pmap(
  list(.x=plot_out$Pred_obs,
       .y=plot_out$Endpoint,
       .z=plot_out$trans,
       .a=plot_out$brks
  ),
  function(.x,.y,.z,.a) 
    .x %>% 
    mutate(WSHED2=WSHED) %>% 
    group_by(WSHED2) %>% 
    nest() %>% 
    ungroup() %>% 
    mutate(plt = map2(data,.y,
                      function(x,y) {
                        ggplot(
                          x,
                          aes(x = Year_n,
                              colour = SiteCode,
                              group = SiteCode
                          )
                        ) +
                          geom_linerange(aes(ymin=pred_low,ymax=pred_high),
                                         colour='grey',
                                         alpha=0.40,
                                         linewidth = 5) +
                          geom_point(aes(y=value),size=3) +
                          scale_colour_discrete() +
                          theme_bw() +
                          ylab(paste0(y)) +
                          xlab("")+
                          facet_wrap(~SiteCode) +
                          labs(
                            title = paste0(x$WSHED[[1]]," | ",.y),
                            subtitle = "Points represent observed values, grey bars represent 90 percent confidence interval from GLMM model"
                          )+
                          theme(
                            axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1),
                            legend.position = "none"
                          )
                      })) %>% 
    pull(plt) %>% 
    map(~.x + scale_y_continuous(breaks = .a(), expand = c(0,0)))%>% 
    map(~.x + coord_transform(y = .z))
)

saveRDS(plot_out,"data/rawdata_plots.rds")

plot_out <- plot_out %>% 
  select(-raw_plots)

# Model fit plot ----------------------------------------------------------

plot_out$model_fit <- pmap(
  list(.x=plot_out$Pred_obs,
       .y=plot_out$Endpoint,
       .z=plot_out$trans,
       .a=plot_out$brks
  ),
  function(.x,.y,.z,.a) {
    ggplot(
      .x,
      aes(x = pred.fit,
          y = value
      )
    ) +
      geom_point(colour='grey25') +
      geom_abline(slope = 1, intercept = 0, linewidth = 0.25,colour = "black")+
      theme_bw() +
      ylab("Observed") +
      xlab("Model Fitted") +
      scale_x_continuous(limits = range(.a(5)(c(.x$pred.fit,.x$value))),breaks = .a()) +
      scale_y_continuous(limits = range(.a(5)(c(.x$pred.fit,.x$value))),breaks = .a()) +
      coord_transform(
        x = .z,
        y = .z
      ) + 
      labs(
        title = paste0(.y),
        subtitle = paste0("Pseudo R2 = ", scales::label_number(0.001)(cor(.x$pred.fit,.x$value,method="spearman")^2))
      )
  }
)


# Model Residuals --------------------------------------------------------------

plot_out$fit1 <- map2(mod_data_fin$sim_resid_refitF,
                      mod_data_fin$Endpoint,
                      function(x,y){
                        ggplotify::as.ggplot(function() plot(x,title = paste("Model Fits -",y)))
                      }
)


# Temporal Plots ----------------------------------------------------------

plot_out$temporal_plot1 <- pmap(
  list(.x=plot_out$Pred_current,
       .y=plot_out$Endpoint,
       .z=plot_out$trans,
       .a=plot_out$brks
  ),
  function(.x,.y,.z,.a) {
    ggplot(
      .x,
      aes(x = SiteCode,
          y = mean_slope,
          colour = WSHED,
          alpha = case_when(
            grepl("Sig",linear_trend) ~ 1,
            grepl("Trend",linear_trend) ~ 0.92,
            T ~ 0.8
          )
      )
    ) +
      geom_pointrange(aes(ymin = lower_slope_sig,
                          ymax = upper_slope_sig,
      ),linewidth = 0.25) +
      geom_pointrange(aes(ymin = lower_slope_trend,
                          ymax = upper_slope_trend,
      ),linewidth = 1.5) +
      geom_point(aes(x = SiteCode,
                     y = mean_slope),
                 shape=21,size=2.5,inherit.aes = F)+
      geom_hline(yintercept = 0) +
      scale_colour_discrete() +
      facet_wrap(~WSHED, scales = "free_x") +
      theme_bw() +
      ylab(.y) +
      theme(
        axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1),
        legend.position = "none"
      ) +
      labs(title = paste0(.y),
           subtitle  = "Mean +/- 50 and 90 percent confidence intervals of the slope coefficients from GLMM models",
           caption = "Areas with confidence intervals that don't overlap 0 are considered to increase or decrease significantly."
      ) 
  }
)

plot_out$temporal_plot1_1 <- pmap(
  list(.x=plot_out$Pred_current,
       .y=plot_out$Endpoint,
       .z=plot_out$trans,
       .a=plot_out$brks
  ),
  function(.x,.y,.z,.a) {
    ggplot(
      .x,
      aes(x = SiteCode,
          y = mean_slope,
          alpha = case_when(
            grepl("Sig",linear_trend) ~ 1,
            grepl("Trend",linear_trend) ~ 0.92,
            T ~ 0.8
          )
      )
    ) +
      geom_linerange(aes(ymin = lower_slope_sig,
                         ymax = upper_slope_sig,
                         colour = `UCA Class`
      ),linewidth = 0.25) +
      geom_linerange(aes(ymin = lower_slope_trend,
                         ymax = upper_slope_trend,
                         colour = `UCA Class`
      ),linewidth = 1.5) +
      geom_point(aes(x = SiteCode,
                     y = mean_slope,
                     shape = `Stress Class`,
                     fill = `UCA Class`,
                     alpha = case_when(
                       grepl("Sig",linear_trend) ~ 1,
                       grepl("Trend",linear_trend) ~ 0.92,
                       T ~ 0.8
                     )),
                 colour = "black",
                 size=2.5,inherit.aes = F)+
      geom_hline(yintercept = 0) +
      scale_shape_manual(values=c(21,22,23)) +
      scale_colour_manual(values = cols_list) +
      scale_fill_manual(values = cols_list) +
      facet_wrap(~WSHED, scales = "free_x") +
      theme_bw() +
      ylab(.y) +
      theme(
        axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1),
        legend.position = "top"
      ) +
      labs(title = paste0(.y),
           subtitle  = "Mean +/- 50 and 90 percent confidence intervals of the slope coefficients from GLMM models",
           caption = "Areas with confidence intervals that don't overlap 0 are considered to increase or decrease significantly."
      ) +
      guides(alpha="none")
  }
)


# Temporal plots by station -----------------------------------------------

plot_out_stntemp <- plot_out %>% 
  select(TRCA,Endpoint,Endpoint_group,fam,trans,brks,Pred_current) %>% 
  unnest(Pred_current) %>% 
  group_by(TRCA) %>% 
  nest() %>% 
  mutate(
    temporal_plot1_st = map(data,
                            ~{
                              ggplot(
                                .x 
                                ,
                                aes(x = Endpoint,
                                    y = mean_slope#,
                                    # alpha = case_when(
                                    #   grepl("Sig",linear_trend) ~ 1,
                                    #   grepl("Trend",linear_trend) ~ 0.92,
                                    #   T ~ 0.8
                                    # )
                                )
                              ) +
                                geom_linerange(aes(ymin = lower_slope_sig,
                                                   ymax = upper_slope_sig,
                                                   colour = Endpoint_group
                                ),linewidth = 0.25) +
                                geom_linerange(aes(ymin = lower_slope_trend,
                                                   ymax = upper_slope_trend,
                                                   colour = Endpoint_group
                                ),linewidth = 1.5) +
                                geom_point(aes(x = Endpoint,
                                               y = mean_slope,
                                               fill = Endpoint_group#,
                                               # alpha = case_when(
                                               #   grepl("Sig",linear_trend) ~ 1,
                                               #   grepl("Trend",linear_trend) ~ 0.92,
                                               #   T ~ 0.8
                                               # )
                                ),
                                colour = "black",
                                shape=21,
                                size=2.5,inherit.aes = F)+
                                geom_hline(yintercept = 0) +
                                facet_wrap(~WSHED+SiteCode, scales = "free",ncol=4) +
                                theme_bw() +
                                ylab("") +
                                xlab("")+
                                theme(
                                  axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1),
                                  legend.position = "none"
                                ) +
                                labs(#title = paste0(.y),
                                  subtitle  = "Mean +/- 50 and 90 percent confidence intervals of the slope coefficients from GLMM models",
                                  caption = "Areas with confidence intervals that don't overlap 0 are considered to increase or decrease significantly."
                                ) +
                                guides(alpha="none")
                              
                            })
  )

saveRDS(plot_out_stntemp,"data/plot_out_stntemp.rds")

# Temporal Plots - High Level ---------------------------------------------

plot_out$temporal_plot2 <- pmap(
  list(.x=plot_out$summary_terms_temporal,
       .y=plot_out$Endpoint,
       .z=plot_out$trans,
       .a=plot_out$brks
  ),
  function(.x,.y,.z,.a) {
    .x %>% 
      mutate(
        gp = case_when(
          group_var == "`1`" ~ "Overall",
          #group_var == "`Habitat Stability Class`" ~ `Habitat Stability Class`,
          group_var == "`Stress Class`" ~ `Stress Class`,
          group_var == "`UCA Class`" ~ `UCA Class`,
        ),
        gp = factor(gp,levels=c("Overall","Small","Intermediate","Large","High","Medium","Low"))
      ) %>% 
      filter(group_var == "`UCA Class`:`Stress Class`") %>% #:`Habitat Stability Class`
      mutate(`Stress Class` = paste0("Stress - ",`Stress Class`)) %>%
      #mutate(`Habitat Stability Class` = paste0("Habitat Stability - ",`Habitat Stability Class`)) %>%
      mutate(`Stress Class` = factor(`Stress Class`,levels = paste0("Stress - ",c("High","Medium","Low")))) %>% 
      #mutate(`Habitat Stability Class` = factor(`Habitat Stability Class`,levels = paste0("Habitat Stability - ",c("High","Low")))) %>% 
      mutate(gp=case_when(
        #is.na(`Stress Class`) ~ `Habitat Stability Class`,
        T ~ `Stress Class`)) %>% 
      ggplot(
        aes(x = `UCA Class`,
            y = mean_slope,
            colour = `Stress Class`, # = `Habitat Stability Class`,
            group = `Stress Class`,# = `Habitat Stability Class`,
            # alpha = case_when(
            #   grepl("Sig",linear_trend) ~ 1,
            #   grepl("Trend",linear_trend) ~ 0.92,
            #   T ~ 0.8
            # )
        )
      ) +
      geom_pointrange(aes(ymin = lower_slope_sig,
                          ymax = upper_slope_sig,
      ),linewidth = 0.25,position=position_dodge(width=0.5)) +
      geom_pointrange(aes(ymin = lower_slope_trend,
                          ymax = upper_slope_trend,
      ),linewidth = 1.5,position=position_dodge(width=0.5)) +
      geom_point(aes(x = `UCA Class`,
                     y = mean_slope,
                     group = `Stress Class`), #`Habitat Stability Class`
                 shape=21,size=2.5,inherit.aes = F,position=position_dodge(width=0.5))+
      geom_hline(yintercept = 0) +
      scale_colour_discrete() +
      facet_grid(~`Stress Class`) +
      theme_bw() +
      ylab(.y) +
      theme(
        legend.position = "bottom"
      ) +
      labs(title = paste0(.y, " by Stress Class"), #and Habitat Stability 
           subtitle  = "Mean +/- 50 and 90 percent confidence intervals of the slope coefficients from GLMM models",
           caption = "Coefficients with confidence intervals that don't overlap 0 are considered to increase or decrease significantly."
      ) +
      guides(alpha="none")
    
  }
)

# Version 2

temp_summary_plot1 <- plot_out %>%
  select(TRCA,
         summary_terms_temporal,
         Endpoint,
         trans,
         brks) %>% 
  unnest(summary_terms_temporal) %>% 
  filter(group_var %in% c("`UCA Class`:`Stress Class`")) %>% #"`Habitat Stability Class`",
  group_by(TRCA,Endpoint) %>% 
  nest() %>% 
  mutate(plot = map2(data,Endpoint,
                     ~{
                       .x  %>% 
                         ggplot(
                           aes(x = `UCA Class`,
                               y = mean_slope,
                               colour = `Stress Class`#,
                               # alpha = case_when(
                               #   grepl("Sig",linear_trend) ~ 1,
                               #   grepl("Trend",linear_trend) ~ 0.92,
                               #   T ~ 0.8
                               # )
                           )
                         ) +
                         geom_pointrange(aes(ymin = lower_slope_sig,
                                             ymax = upper_slope_sig,
                         ),linewidth = 0.25,
                         position=position_dodge(width=0.25)) +
                         geom_pointrange(aes(ymin = lower_slope_trend,
                                             ymax = upper_slope_trend,
                         ),linewidth = 1.5,
                         position=position_dodge(width=0.25)) +
                         geom_point(aes(x = `UCA Class`,
                                        y = mean_slope,
                                        group = `Stress Class`),
                                    shape=21,size=2.5,inherit.aes = F,
                                    position=position_dodge(width=0.25)) +
                         geom_hline(yintercept = 0) +
                         scale_colour_discrete() +
                         theme_bw() +
                         ylab("")+
                         xlab("")+
                         theme(
                           legend.position = "bottom"
                         ) +
                         labs(title = .y)+
                         guides(alpha="none")
                     }))

temp_summary_plot1 <- temp_summary_plot1 %>% 
  ungroup() %>% 
  group_by(TRCA) %>% 
  nest() %>% 
  mutate(plot = map(data,~cowplot::plot_grid(plotlist=.x$plot,ncol=1,align ="hv",axis ="tblr"))) %>% 
  mutate(plot = map(plot,~cowplot::plot_grid( 
    cowplot::ggdraw()+cowplot::draw_label(
      "Overall Trends",
      fontface = 'bold',
      size = 14,
      x = 0,
      hjust = 0
    ),
    cowplot::ggdraw()+cowplot::draw_label(
      "Mean +/- 50 and 90 percent confidence intervals of the slope coefficients from GLMM models",
      fontface = 'plain',
      size = 12,
      x = 0,
      hjust = 0
    ),
    cowplot::ggdraw()+cowplot::draw_label(
      "Coefficients with confidence intervals that don't overlap 0 are considered to increase or decrease significantly",
      fontface = 'plain',
      size = 10,
      x = 0,
      hjust = 0
    ),
    .x,
    ncol=1,
    rel_heights = c(0.01,0.01,0.01,1))))

saveRDS(temp_summary_plot1,"data/temp_summary_plot1.rds")


# Temporal Plots - Summary ------------------------------------------------

temp_summary_plot <- plot_out %>%
  select(TRCA,
         summary_terms_temporal,
         Endpoint,
         trans,
         brks) %>% 
  unnest(summary_terms_temporal) %>% 
  filter(group_var %in% c("`Stress Class`","`UCA Class`","`1`")) %>% #"`Habitat Stability Class`",
  mutate(
    gp = case_when(
      group_var == "`1`" ~ "Overall",
      #group_var == "`Habitat Stability Class`" ~ `Habitat Stability Class`,
      group_var == "`Stress Class`" ~ `Stress Class`,
      group_var == "`UCA Class`" ~ `UCA Class`,
    ),
    group_var = case_when(
      group_var == "`1`" ~ "Overall",
      #group_var == "`Habitat Stability Class`" ~ "Habitat Stability Class",
      group_var == "`Stress Class`" ~ "Stress Class",
      group_var == "`UCA Class`" ~ "UCA Class",
      T ~ group_var
    )
  ) %>% 
  mutate(group_var = factor(group_var, levels = rev(c("Stress Class","UCA Class","Overall")))) %>% #"Habitat Stability Class",
  group_by(TRCA,Endpoint) %>% 
  nest() %>% 
  mutate(plot = map2(data,Endpoint,
                     ~{
                       .x  %>% 
                         mutate(gp = factor(gp,levels=c("Overall","Small","Intermediate","Large","High","Medium","Low"))) %>% 
                         #filter(!is.na(gp)) %>% 
                         ggplot(
                           aes(x = gp,
                               y = mean_slope,
                               colour = group_var#,
                               # alpha = case_when(
                               #   grepl("Sig",linear_trend) ~ 1,
                               #   grepl("Trend",linear_trend) ~ 0.92,
                               #   T ~ 0.8
                               # )
                           )
                         ) +
                         geom_pointrange(aes(ymin = lower_slope_sig,
                                             ymax = upper_slope_sig,
                         ),linewidth = 0.25) +
                         geom_pointrange(aes(ymin = lower_slope_trend,
                                             ymax = upper_slope_trend,
                         ),linewidth = 1.5) +
                         geom_point(aes(x = gp,
                                        y = mean_slope),
                                    shape=21,size=2.5,inherit.aes = F) +
                         geom_hline(yintercept = 0) +
                         scale_colour_discrete() +
                         facet_grid(~group_var,scales="free", drop=TRUE) +
                         theme_bw() +
                         ylab("")+
                         xlab("")+
                         theme(
                           legend.position = "none"
                         ) +
                         labs(title = .y)+
                         guides(alpha="none")
                     }))

temp_summary_plot <- temp_summary_plot %>% 
  ungroup() %>% 
  group_by(TRCA) %>% 
  nest() %>% 
  mutate(plot = map(data,~cowplot::plot_grid(plotlist=.x$plot,ncol=1,align ="hv",axis ="tblr"))) %>% 
  mutate(plot = map(plot,~cowplot::plot_grid( 
    cowplot::ggdraw()+cowplot::draw_label(
      "Overall Trends",
      fontface = 'bold',
      size = 14,
      x = 0,
      hjust = 0
    ),
    cowplot::ggdraw()+cowplot::draw_label(
      "Mean +/- 50 and 90 percent confidence intervals of the slope coefficients from GLMM models",
      fontface = 'plain',
      size = 12,
      x = 0,
      hjust = 0
    ),
    cowplot::ggdraw()+cowplot::draw_label(
      "Coefficients with confidence intervals that don't overlap 0 are considered to increase or decrease significantly",
      fontface = 'plain',
      size = 10,
      x = 0,
      hjust = 0
    ),
    .x,
    ncol=1,
    rel_heights = c(0.01,0.01,0.01,1))))

saveRDS(temp_summary_plot,"data/temp_summary_plot.rds")

# Temporal Maps -----------------------------------------------------------

plot_out$temporal_map1 <- map2(
  plot_out$Pred_current,
  plot_out$Endpoint,
  ~basemap_ggplot(sf::st_bbox(TRCA_subwats),map_service="carto",map_type="light",verbose=F) +
    theme_minimal()+
    theme(legend.position = "bottom",
          axis.text = element_blank(),
          axis.ticks = element_blank(),
          legend.key.width = unit(dev.size()[1] / 15, "inches"),
          legend.key.height = unit(dev.size()[2] / 50, "inches"),
          plot.margin=unit(c(0,0,0,0),"cm")) +
    ylab("")+
    xlab("")+
    geom_sf(data=sf::st_transform(TRCA_subwats,crs="EPSG: 3857"),alpha=0.75)+
    geom_sf(data=sf::st_transform(TRCA_wats,crs="EPSG: 3857"),alpha=0,linewidth=1.25)+
    geom_sf(data=sf::st_transform(strm_lns,crs="EPSG: 3857"),alpha=0.25,colour="darkblue")+
    geom_sf(data=sf::st_transform(.x,crs="EPSG: 3857"),
            aes(colour=linear_trend),
            size=3,
            linewidth=0.05) +
    geom_sf(data=sf::st_transform(.x,crs="EPSG: 3857"),
            size=3,
            shape=21) +
    scale_colour_manual(values = c(
      `Increasing Significant`="darkolivegreen4",
      `Increasing Trend`="darkolivegreen2",
      `No Change`="black",
      `Decreasing Trend`="skyblue2",
      `Decreasing Significant`="skyblue4"
    )) +
    labs(
      title = .y,
      subtitle  = "Mean slope coefficients from GLMM models",
      colour = .y
    ) 
)



# Spatial Plots -----------------------------------------------------------


plot_out$spatial_plot1 <- pmap(
  list(.x=plot_out$Pred_current,
       .y=plot_out$Endpoint,
       .z=plot_out$trans,
       .a=plot_out$brks
  ),
  function(.x,.y,.z,.a) {
    ggplot(
      .x,
      aes(x = SiteCode,
          y = pred.fit,
          colour = WSHED
      )
    ) +
      geom_pointrange(aes(ymin = pred_low,
                          ymax = pred_high,
      ),linewidth = 0.25) +
      geom_pointrange(aes(ymin = pred_low_trend,
                          ymax = pred_high_trend,
      ),linewidth = 1.5) +
      geom_point(aes(x = SiteCode,
                     y = pred.fit),shape=21,size=2.5,inherit.aes = F)+
      
      scale_colour_discrete() +
      facet_wrap(~WSHED, scales = "free_x") +
      theme_bw() +
      ylab(.y) +
      theme(
        axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1),
        legend.position = "none"
      ) +
      labs(title = paste0(.y),
           subtitle = "Mean +/- 50 and 90 percent confidence intervals of predicted 2024 mean values from GLMM models",
           caption  = "Areas with confidence intervals that don't overlap 0 are considered to differ significantly"
      ) +
      scale_y_continuous(breaks = .a()) +
      coord_transform(
        y = .z
      )
  }
)

plot_out$spatial_plot1_1 <- pmap(
  list(.x=plot_out$Pred_current,
       .y=plot_out$Endpoint,
       .z=plot_out$trans,
       .a=plot_out$brks
  ),
  function(.x,.y,.z,.a) {
    ggplot(
      .x,
      aes(x = SiteCode,
          y = pred.fit,
          shape = `Stress Class`,
          colour = `UCA Class`
      )
    ) +
      geom_linerange(aes(ymin = pred_low,
                         ymax = pred_high,
                         colour = `UCA Class`
      ),linewidth = 0.25) +
      geom_linerange(aes(ymin = pred_low_trend,
                         ymax = pred_high_trend,
                         colour = `UCA Class`
      ),linewidth = 1.5) +
      geom_point(aes(x = SiteCode,
                     y = pred.fit,
                     shape = `Stress Class`,
                     fill = `UCA Class`),size=2.5,inherit.aes = F)+
      scale_shape_manual(values=c(21,22,23))+
      scale_colour_manual(values = cols_list) +
      scale_fill_manual(values = cols_list) +
      facet_wrap(~WSHED, scales = "free_x") +
      theme_bw() +
      ylab(.y) +
      theme(
        axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1),
        legend.position = "top"
      ) +
      labs(title = paste0(.y),
           subtitle = "Mean +/- 50 and 90 percent confidence intervals of predicted 2024 mean values from GLMM models",
           caption  = "Areas with confidence intervals that don't overlap 0 are considered to differ significantly"
      ) +
      scale_y_continuous(breaks = .a()) +
      coord_transform(
        y = .z
      )
  }
)


# Spatial Plot - High Level -----------------------------------------------


plot_out$spatial_plot2 <- pmap(
  list(.x=plot_out$summary_terms_spatial,
       .y=plot_out$Endpoint,
       .z=plot_out$trans,
       .a=plot_out$brks
  ),
  function(.x,.y,.z,.a) {
    .x %>% 
      mutate(
        gp = case_when(
          group_var == "`1`" ~ "Overall",
          #group_var == "`Habitat Stability Class`" ~ `Habitat Stability Class`,
          group_var == "`Stress Class`" ~ `Stress Class`,
          group_var == "`UCA Class`" ~ `UCA Class`,
        ),
        gp = factor(gp,levels=c("Overall","Small","Intermediate","Large","High","Medium","Low"))
      ) %>% 
      filter(group_var == "`UCA Class`:`Stress Class`") %>% #:`Habitat Stability Class`
      mutate(`Stress Class` = paste0("Stress - ",`Stress Class`)) %>%
      #mutate(`Habitat Stability Class` = paste0("Habitat Stability - ",`Habitat Stability Class`)) %>%
      mutate(`Stress Class` = factor(`Stress Class`,levels = paste0("Stress - ",c("High","Medium","Low")))) %>% 
      #mutate(`Habitat Stability Class` = factor(`Habitat Stability Class`,levels = paste0("Habitat Stability - ",c("High","Low")))) %>% 
      mutate(gp=case_when(
        #is.na(`Stress Class`) ~ `Habitat Stability Class`,
        T ~ `Stress Class`)) %>% 
      ggplot(
        aes(x = `UCA Class`,
            y = value,
            colour = `Stress Class`, #`Habitat Stability Class`
            group = `Stress Class` #`Habitat Stability Class`
        )
      ) +
      geom_pointrange(aes(ymin = pred_low,
                          ymax = pred_high,
      ),linewidth = 0.25,position=position_dodge(width=0.5)) +
      geom_pointrange(aes(ymin = pred_low_trend,
                          ymax = pred_high_trend,
      ),linewidth = 1.5,position=position_dodge(width=0.5)) +
      geom_point(aes(x = `UCA Class`,
                     y = value,
                     group = `Stress Class`), #`Habitat Stability Class`
                 shape=21,size=2.5,inherit.aes = F,position=position_dodge(width=0.5))+
      
      scale_colour_discrete() +
      #facet_grid(~`Stress Class`) +
      theme_bw() +
      ylab(.y) +
      theme(
        legend.position = "bottom"
      ) +
      labs(title = paste0(.y),
           subtitle = "Mean +/- 50 and 90 percent confidence intervals of predicted 2024 mean values from GLMM models",
           caption  = "Coefficients with confidence intervals that don't overlap are considered to differ significantly"
      ) +
      scale_y_continuous(breaks = .a()) +
      coord_transform(
        y = .z
      )
  }
)

#version 2
spat_summary_plot1 <- plot_out %>%
  select(TRCA,
         summary_terms_spatial,
         Endpoint,
         trans,
         brks) %>% 
  unnest(summary_terms_spatial) %>% 
  filter(group_var %in% c("`UCA Class`:`Stress Class`")) %>% #"`Habitat Stability Class`",
  group_by(TRCA,Endpoint) %>% 
  nest() %>% 
  mutate(plot = map2(data,Endpoint,
                     ~{
                       .x  %>% 
                         ggplot(
                           aes(x = `UCA Class`,
                               y = value,
                               colour = `Stress Class`
                           )
                         ) +
                         geom_pointrange(
                           aes(ymin = pred_low,
                               ymax = pred_high,
                           ),
                           linewidth = 0.25,
                           position=position_dodge(width=0.25)) +
                         geom_pointrange(
                           aes(ymin = pred_low_trend,
                               ymax = pred_high_trend,
                           ),
                           linewidth = 1.5,
                           position=position_dodge(width=0.25)) +
                         geom_point(
                           aes(x = `UCA Class`,
                               y = value,
                               group = `Stress Class`),
                           shape=21,
                           size=2.5,
                           inherit.aes = F,
                           position=position_dodge(width=0.25)) +
                         scale_colour_discrete() +
                         theme_bw() +
                         ylab("")+
                         xlab("")+
                         theme(
                           legend.position = "bottom"
                         ) +
                         labs(title=.y)+
                         guides(alpha="none") +
                         scale_y_continuous(breaks = .x$brks[[1]]()) +
                         coord_transform(
                           y = .x$trans[[1]]
                         )
                     }))

spat_summary_plot1 <- spat_summary_plot1 %>% 
  ungroup() %>% 
  group_by(TRCA) %>% 
  nest() %>% 
  mutate(plot = map(data,~cowplot::plot_grid(plotlist=.x$plot,ncol=1,align ="hv",axis ="tblr") )) %>% 
  mutate(plot = map(plot,~cowplot::plot_grid(
    cowplot::ggdraw()+cowplot::draw_label(
      "Overall Means",
      fontface = 'bold',
      size = 14,
      x = 0,
      hjust = 0
    ),
    cowplot::ggdraw()+cowplot::draw_label(
      "Mean +/- 50 and 90 percent confidence intervals of predicted 2024 mean values from GLMM models",
      fontface = 'plain',
      size = 12,
      x = 0,
      hjust = 0
    ),
    cowplot::ggdraw()+cowplot::draw_label(
      "Coefficients with confidence intervals that don't overlap are considered to differ significantly",
      fontface = 'plain',
      size = 10,
      x = 0,
      hjust = 0
    ),
    .x,
    ncol=1,
    rel_heights = c(0.01,0.01,0.01,1))))

saveRDS(spat_summary_plot1,"data/spat_summary_plot1.rds")


# Spatial Plots - Summary ------------------------------------------------

spat_summary_plot <- plot_out %>%
  select(TRCA,
         summary_terms_spatial,
         Endpoint,
         trans,
         brks) %>% 
  unnest(summary_terms_spatial) %>% 
  filter(group_var %in% c("`Stress Class`","`UCA Class`","`1`")) %>% #"`Habitat Stability Class`",
  mutate(
    gp = case_when(
      group_var == "`1`" ~ "Overall",
      #group_var == "`Habitat Stability Class`" ~ `Habitat Stability Class`,
      group_var == "`Stress Class`" ~ `Stress Class`,
      group_var == "`UCA Class`" ~ `UCA Class`,
    ),
    group_var = case_when(
      group_var == "`1`" ~ "Overall",
      #group_var == "`Habitat Stability Class`" ~ "Habitat Stability Class",
      group_var == "`Stress Class`" ~ "Stress Class",
      group_var == "`UCA Class`" ~ "UCA Class",
      T ~ group_var
    )
  ) %>% 
  mutate(group_var = factor(group_var, levels = rev(c("Stress Class","UCA Class","Overall")))) %>% #"Habitat Stability Class",
  group_by(TRCA,Endpoint) %>% 
  nest() %>% 
  mutate(plot = map2(data,Endpoint,
                     ~{
                       .x  %>% 
                         mutate(gp = factor(gp,levels=c("Overall","Small","Intermediate","Large","High","Medium","Low"))) %>% 
                         ggplot(
                           aes(x = gp,
                               y = value,
                               colour = group_var
                           )
                         ) +
                         geom_pointrange(aes(ymin = pred_low,
                                             ymax = pred_high,
                         ),linewidth = 0.25) +
                         geom_pointrange(aes(ymin = pred_low_trend,
                                             ymax = pred_high_trend,
                         ),linewidth = 1.5) +
                         geom_point(aes(x = gp,
                                        y = value),
                                    shape=21,size=2.5,inherit.aes = F) +
                         scale_colour_discrete() +
                         facet_grid(~group_var,scales="free", drop=TRUE) +
                         theme_bw() +
                         ylab("")+
                         xlab("")+
                         theme(
                           legend.position = "none"
                         ) +
                         labs(title=.y)+
                         guides(alpha="none") +
                         scale_y_continuous(breaks = .x$brks[[1]]()) +
                         coord_transform(
                           y = .x$trans[[1]]
                         )
                     }))

spat_summary_plot <- spat_summary_plot %>% 
  ungroup() %>% 
  group_by(TRCA) %>% 
  nest() %>% 
  mutate(plot = map(data,~cowplot::plot_grid(plotlist=.x$plot,ncol=1,align ="hv",axis ="tblr") )) %>% 
  mutate(plot = map(plot,~cowplot::plot_grid(
    cowplot::ggdraw()+cowplot::draw_label(
      "Overall Means",
      fontface = 'bold',
      size = 14,
      x = 0,
      hjust = 0
    ),
    cowplot::ggdraw()+cowplot::draw_label(
      "Mean +/- 50 and 90 percent confidence intervals of predicted 2024 mean values from GLMM models",
      fontface = 'plain',
      size = 12,
      x = 0,
      hjust = 0
    ),
    cowplot::ggdraw()+cowplot::draw_label(
      "Coefficients with confidence intervals that don't overlap are considered to differ significantly",
      fontface = 'plain',
      size = 10,
      x = 0,
      hjust = 0
    ),
    .x,
    ncol=1,
    rel_heights = c(0.01,0.01,0.01,1))))

saveRDS(spat_summary_plot,"data/spat_summary_plot.rds")

# Spatial Maps -----------------------------------------------------------

plot_out$spatial_map1 <- pmap(
  list(.x=plot_out$Pred_current,
       .y=plot_out$Endpoint,
       .z=plot_out$trans,
       .a=plot_out$brks
  ),
  function(.x,.y,.z,.a) {
    basemap_ggplot(sf::st_bbox(TRCA_subwats),map_service="carto",map_type="light",verbose=F) +
      theme_minimal()+
      theme(legend.position = "bottom",
            axis.text = element_blank(),
            axis.ticks = element_blank(),
            legend.key.width = unit(dev.size()[1] / 15, "inches"),
            legend.key.height = unit(dev.size()[2] / 50, "inches"),
            plot.margin=unit(c(0,0,0,0),"cm")) +
      ylab("")+
      xlab("")+
      geom_sf(data=sf::st_transform(TRCA_subwats,crs="EPSG: 3857"),alpha=0.75)+
      geom_sf(data=sf::st_transform(TRCA_wats,crs="EPSG: 3857"),alpha=0,linewidth=1.25)+
      geom_sf(data=sf::st_transform(strm_lns,crs="EPSG: 3857"),alpha=0.25,colour="darkblue")+
      geom_sf(data=sf::st_transform(.x,crs="EPSG: 3857"),
              aes(colour=pred.fit),
              size=3,
              linewidth=0.1) +
      scale_colour_continuous(palette = "Spectral",type="div",trans = .z,breaks=.a()) +
      geom_sf(data=sf::st_transform(.x,crs="EPSG: 3857"),
              shape=21,
              size=3) +
      labs(
        title = .y,
        subtitle = "Predicted 2024 mean values from GLMM models",
        colour = .y
      ) 
    
  }
)


saveRDS(plot_out,"data/results_plots.rds")



# Temporal Summary Tables ----------------------------------------------------------

tbl_out <- plot_out %>% 
  select(Endpoint_group,Endpoint,Pred_current) %>% 
  unnest(Pred_current) %>% 
  mutate(`UCA Class`=factor(`UCA Class`,levels = c("Small","Intermediate","Large"))) %>% 
  mutate(`Stress Class`=factor(`Stress Class`,levels = c("Low","Medium","High"))) %>% 
  #mutate(`Habitat Stability Class`=factor(`Habitat Stability Class`,levels = c("High","Low"))) %>% 
  arrange(`Stress Class`,`UCA Class`) %>% #,`Habitat Stability Class`
  mutate(value = paste0(mean_slope," (",lower_slope_sig," - ",upper_slope_sig,")")) %>% 
  select(Endpoint_group,Endpoint,WSHED,SiteCode,value,linear_trend,
         `Stress Class`,`UCA Class`#,`Habitat Stability Class`
  )

ep_grp <- tbl_out %>% 
  select(Endpoint_group,Endpoint) %>% 
  distinct()%>% 
  group_by(Endpoint_group) %>% 
  group_split() %>% 
  map(~setNames(list(.x$Endpoint),.x$Endpoint_group[[1]])) %>% 
  unlist(recursive = F)

site_grp <- tbl_out %>% 
  select(WSHED,SiteCode,`Stress Class`,`UCA Class`) %>%  #,`Habitat Stability Class`
  distinct() %>% 
  group_by(WSHED) %>% 
  group_split() %>% 
  map(~setNames(list(.x %>% select(-WSHED)),.x$WSHED[[1]])) %>% 
  unlist(recursive = F)

shade_tbl <- tbl_out %>% 
  select(-value,-Endpoint_group,-WSHED) %>% 
  pivot_wider(names_from = Endpoint,values_from = linear_trend) %>% 
  mutate(SiteCode = as.character(SiteCode)) %>% 
  mutate(across(everything(),~case_when(is.na(.x) ~ "No Change", T ~ .x)))

writexl::write_xlsx(shade_tbl,"data/changetbl1.xlsx")

shade_tbl_sym <- shade_tbl %>% 
  mutate(across(everything(),
                ~case_when(
                  .x == "Increasing Significant" ~ "++",
                  .x == "Increasing Trend" ~ "+",
                  .x == "No Change" ~ "0",
                  .x == "Decreasing Trend" ~ "-",
                  .x == "Decreasing Significant" ~ "--",
                  is.na(.x) ~ "0",
                  T ~ .x
                )))

num_tbl <- tbl_out %>% 
  select(-linear_trend,-Endpoint_group,-WSHED) %>% 
  pivot_wider(names_from = Endpoint,values_from = value) %>% 
  mutate(SiteCode = as.character(SiteCode))

col_list <- list(`Increasing Significant` = "darkolivegreen4",
                 `Increasing Trend` = "darkolivegreen2",
                 `No Change` = "grey90",
                 `Decreasing Trend` = "skyblue2",
                 `Decreasing Significant` = "skyblue4"
)


col_list_leg <- list(`+ + = Increasing Significant` = "darkolivegreen4",
                     `+ = Increasing Trend` = "darkolivegreen2",
                     `0 = No Change` = "darkgrey",
                     `- = Decreasing Trend` = "skyblue2",
                     `- - = Decreasing Significant` = "skyblue4"
)

col_list_leg <- map2(names(col_list_leg),col_list,
                     ~paste0("<span style='color:",col2hex(.y),";'>",.x,"</span>"))

col_list_leg <- paste0(col_list_leg, collapse = " | ")


gt_list <- list()

for (v in names(site_grp)){
  
  shade_tbl_sub <- shade_tbl %>% 
    filter(SiteCode %in% site_grp[[v]]$SiteCode) %>%
    select(-`Stress Class`)
  shade_tbl_sym_sub <- shade_tbl_sym %>%
    filter(SiteCode %in% site_grp[[v]]$SiteCode) %>%
    select(-`Stress Class`)
  
  gt_tbl <- gt::gt(shade_tbl_sym_sub  ,
                   rowname_col = "SiteCode") %>% 
    gt::tab_header(
      title = v
    ) %>%
    gt::tab_stubhead(label = v) %>% 
    gt::tab_source_note(
      source_note = gt::md(col_list_leg)
    ) %>% 
    gt::tab_style(
      style = gt::cell_text(align = "center", v_align = "middle"),
      locations = gt::cells_body(
        columns = everything(),
        rows = everything()
      )
    ) %>% 
    gt::tab_style(
      style = gt::cell_text(align = "center",
                            v_align = "middle",
                            stretch = "extra-condensed",
                            size = "small"),
      locations = gt::cells_column_labels()
    ) %>% 
    gt::tab_style(
      style = gt::cell_text(align = "center",
                            v_align = "middle",
                            stretch = "extra-condensed",
                            size = "small",
                            transform="capitalize"),
      locations = gt::cells_row_groups()
    )
  
  for (i in names(col_list)){
    for (ii in 1:ncol(shade_tbl_sub)){
      gt_tbl <- gt_tbl %>% 
        gt::tab_style(
          style = gt::cell_fill(color = col_list[[i]]),
          locations = gt::cells_body(
            columns = ii,
            rows = which(shade_tbl_sub[[ii]] == i)
          )
        )
    }
  }
  
  for (i in names(ep_grp)){
    gt_tbl <- gt_tbl %>% gt::tab_spanner(label = i,columns  = as.character(ep_grp[[i]]))
  }
  
  for (i in unique(site_grp[[v]]$`Stress Class`)){
    gt_tbl <- gt_tbl %>% gt::tab_row_group(label = paste0(i," Stess Class"),rows = as.character(site_grp[[v]]$SiteCode[site_grp[[v]]$`Stress Class`==i]))
  }
  
  gt_list[[v]] <- gt_tbl
}

gt_tbl <- gt::gt_group(
  .list = gt_list
)

saveRDS(gt_tbl,"data/summarytbl.rds")

# Trend Summary Table 1 ---------------------------------------------------

tbl_out <- plot_out %>% 
  select(Endpoint_group,Endpoint,summary_terms_temporal) %>% 
  unnest(summary_terms_temporal) %>% 
  filter(group_var %in% c("`1`","`Stress Class`","`UCA Class`")) %>% #,"`Habitat Stability Class`",
  select(Endpoint_group,Endpoint,`Stress Class`,`UCA Class`,`1`,linear_trend) %>% #,`Habitat Stability Class`
  arrange(`Stress Class`,`UCA Class`,`1`) %>% #,`Habitat Stability Class`
  pivot_longer(c(`Stress Class`,`UCA Class`,`1`)) %>% #,`Habitat Stability Class`
  filter(!is.na(value))

levels(tbl_out$value)[levels(tbl_out$value)=="overall"]<-"Global Trend"

ep_grp <- tbl_out %>% 
  select(Endpoint_group,Endpoint) %>% 
  distinct()%>% 
  group_by(Endpoint_group) %>% 
  group_split() %>% 
  map(~setNames(list(.x$Endpoint),.x$Endpoint_group[[1]])) %>% 
  unlist(recursive = F)

site_grp <- tbl_out %>% 
  select(name,value) %>% 
  distinct() %>% 
  group_by(name) %>% 
  group_split() %>% 
  map(~setNames(list(.x$value),.x$name[[1]])) %>% 
  unlist(recursive = F)

shade_tbl <- tbl_out %>% 
  select(-Endpoint_group) %>% 
  pivot_wider(names_from = Endpoint,values_from = linear_trend) %>% 
  arrange(name,value) %>% 
  mutate(across(everything(),~case_when(is.na(.x) ~ "No Change", T ~ .x)))

site_grp <- shade_tbl %>% 
  mutate(value = row_number()) %>% 
  select(name,value) %>% 
  distinct() %>% 
  group_by(name) %>% 
  group_split() %>% 
  map(~setNames(list(.x$value),.x$name[[1]])) %>% 
  unlist(recursive = F)

site_grp<-rev(site_grp)
names(site_grp)[names(site_grp)=="1"] <- " "


shade_tbl <- shade_tbl %>% 
  select(-name) 

shade_tbl_sym <- shade_tbl %>% 
  mutate(across(everything(),
                ~case_when(
                  .x == "Increasing Significant" ~ "++",
                  .x == "Increasing Trend" ~ "+",
                  .x == "No Change" ~ "0",
                  .x == "Decreasing Trend" ~ "-",
                  .x == "Decreasing Significant" ~ "--",
                  is.na(.x) ~ "0",
                  T ~ .x
                )))

gt_tbl1 <- gt::gt(shade_tbl_sym,rowname_col = "value") %>% 
  gt::tab_stubhead(label = "Overall Trends") %>% 
  gt::tab_source_note(
    source_note = gt::md(col_list_leg)
  ) %>% 
  gt::tab_style(
    style = gt::cell_text(align = "center", v_align = "middle"),
    locations = gt::cells_body(
      columns = everything(),
      rows = everything()
    )
  ) %>% 
  gt::tab_style(
    style = gt::cell_text(align = "center",
                          v_align = "middle",
                          stretch = "extra-condensed",
                          size = "small"),
    locations = gt::cells_column_labels()
  ) %>% 
  gt::tab_style(
    style = gt::cell_text(align = "center",
                          v_align = "middle",
                          stretch = "extra-condensed",
                          size = "small",
                          transform="capitalize"),
    locations = gt::cells_row_groups()
  )

for (i in names(col_list)){
  for (ii in 1:ncol(shade_tbl)){
    gt_tbl1 <- gt_tbl1 %>% 
      gt::tab_style(
        style = gt::cell_fill(color = col_list[[i]]),
        locations = gt::cells_body(
          columns = ii,
          rows = which(shade_tbl[[ii]] == i)
        )
      )
  }
}

for (i in names(site_grp)){
  gt_tbl1 <- gt_tbl1 %>% gt::tab_row_group(label = i,rows = site_grp[[i]])
}

for (i in names(ep_grp)){
  gt_tbl1 <- gt_tbl1 %>% gt::tab_spanner(label = i,columns  = as.character(ep_grp[[i]]))
}

# Trend Summary Table 2 ---------------------------------------------------

tbl_out <- plot_out %>% 
  select(Endpoint_group,Endpoint,summary_terms_temporal) %>% 
  unnest(summary_terms_temporal) %>% 
  filter(group_var == "`UCA Class`:`Stress Class`") %>% 
  select(Endpoint_group,Endpoint,`Stress Class`,`UCA Class`,linear_trend) %>% 
  filter(!is.na(`Stress Class`)) %>% 
  mutate(`UCA Class` = paste("UCA",`UCA Class`)) %>% 
  mutate(`UCA Class` = factor(`UCA Class`,levels=c("UCA Small","UCA Intermediate","UCA Large"))) %>% 
  arrange(`UCA Class`)

ep_grp <- tbl_out %>% 
  select(Endpoint_group,Endpoint) %>% 
  distinct()%>% 
  group_by(Endpoint_group) %>% 
  group_split() %>% 
  map(~setNames(list(.x$Endpoint),.x$Endpoint_group[[1]])) %>% 
  unlist(recursive = F)

site_grp <- tbl_out %>% 
  select(`Stress Class`,`UCA Class`) %>% 
  distinct() %>% 
  group_by(`Stress Class`) %>% 
  group_split() %>% 
  map(~setNames(list(.x$`UCA Class`),.x$`Stress Class`[[1]])) %>% 
  unlist(recursive = F)

shade_tbl <- tbl_out %>% 
  select(-Endpoint_group) %>% 
  pivot_wider(names_from = Endpoint,values_from = linear_trend) %>% 
  arrange(`Stress Class`,`UCA Class`) %>% 
  mutate(`UCA Class` = as.character(`UCA Class`)) %>% 
  mutate(across(everything(),~case_when(is.na(.x) ~ "No Change", T ~ .x)))

site_grp <- shade_tbl %>% 
  mutate(`UCA Class` = row_number()) %>% 
  select(`Stress Class`,`UCA Class`) %>% 
  distinct() %>% 
  group_by(`Stress Class`) %>% 
  group_split() %>% 
  map(~setNames(list(.x$`UCA Class`),.x$`Stress Class`[[1]])) %>% 
  unlist(recursive = F)

site_grp <- site_grp[c("Low","Medium","High")]
names(site_grp) <- c(c("Low Stress","Medium Stress","High Stress"))

writexl::write_xlsx(shade_tbl,"data/changetbl2.xlsx")

shade_tbl <- shade_tbl %>% 
  select(-`Stress Class`) 

shade_tbl_sym <- shade_tbl %>% 
  mutate(across(everything(),
                ~case_when(
                  .x == "Increasing Significant" ~ "++",
                  .x == "Increasing Trend" ~ "+",
                  .x == "No Change" ~ "0",
                  .x == "Decreasing Trend" ~ "-",
                  .x == "Decreasing Significant" ~ "--",
                  is.na(.x) ~ "0",
                  T ~ .x
                )))

gt_tbl2 <- gt::gt(shade_tbl_sym,rowname_col = "UCA Class") %>% 
  # gt::tab_header(
  #   title = "Temporal Changes across UCA Class"
  # ) %>% 
  gt::tab_stubhead(label = "UCA Class / Stress Class") %>% 
  gt::tab_source_note(
    source_note = gt::md(col_list_leg)
  ) %>% 
  gt::tab_style(
    style = gt::cell_text(align = "center", v_align = "middle"),
    locations = gt::cells_body(
      columns = everything(),
      rows = everything()
    )
  ) %>% 
  gt::tab_style(
    style = gt::cell_text(align = "center",
                          v_align = "middle",
                          stretch = "extra-condensed",
                          size = "small"),
    locations = gt::cells_column_labels()
  ) %>% 
  gt::tab_style(
    style = gt::cell_text(align = "center",
                          v_align = "middle",
                          stretch = "extra-condensed",
                          size = "small",
                          transform="capitalize"),
    locations = gt::cells_row_groups()
  )

for (i in names(col_list)){
  for (ii in 1:ncol(shade_tbl)){
    gt_tbl2 <- gt_tbl2 %>% 
      gt::tab_style(
        style = gt::cell_fill(color = col_list[[i]]),
        locations = gt::cells_body(
          columns = ii,
          rows = which(shade_tbl[[ii]] == i)
        )
      )
  }
}

for (i in names(site_grp)){
  gt_tbl2 <- gt_tbl2 %>% gt::tab_row_group(label = i,rows = site_grp[[i]])
}

for (i in names(ep_grp)){
  gt_tbl2 <- gt_tbl2 %>% gt::tab_spanner(label = i,columns  = as.character(ep_grp[[i]]))
}


# Trend Summary Table 3 ---------------------------------------------------

# tbl_out <- plot_out %>% 
#   select(Endpoint_group,Endpoint,summary_terms_temporal) %>% 
#   unnest(summary_terms_temporal) %>% 
#   filter(group_var == "`UCA Class`:`Habitat Stability Class`") %>% 
#   select(Endpoint_group,Endpoint,`Habitat Stability Class`,`UCA Class`,linear_trend) %>% 
#   filter(!is.na(`Habitat Stability Class`)) %>% 
#   mutate(`UCA Class` = paste("UCA",`UCA Class`)) %>% 
#   mutate(`UCA Class` = factor(`UCA Class`,levels=c("UCA Small","UCA Intermediate","UCA Large"))) %>% 
#   arrange(`UCA Class`)
# 
# 
# ep_grp <- tbl_out %>% 
#   select(Endpoint_group,Endpoint) %>% 
#   distinct()%>% 
#   group_by(Endpoint_group) %>% 
#   group_split() %>% 
#   map(~setNames(list(.x$Endpoint),.x$Endpoint_group[[1]])) %>% 
#   unlist(recursive = F)
# 
# site_grp <- tbl_out %>% 
#   select(`Habitat Stability Class`,`UCA Class`) %>% 
#   distinct() %>% 
#   group_by(`Habitat Stability Class`) %>% 
#   group_split() %>% 
#   map(~setNames(list(.x$`UCA Class`),.x$`Habitat Stability Class`[[1]])) %>% 
#   unlist(recursive = F)
# 
# shade_tbl <- tbl_out %>% 
#   select(-Endpoint_group) %>% 
#   pivot_wider(names_from = Endpoint,values_from = linear_trend) %>% 
#   arrange(`Habitat Stability Class`,`UCA Class`) %>% 
#   mutate(`UCA Class` = as.character(`UCA Class`)) %>% 
#   mutate(across(everything(),~case_when(is.na(.x) ~ "No Change", T ~ .x)))
# 
# site_grp <- shade_tbl %>% 
#   mutate(`UCA Class` = row_number()) %>% 
#   select(`Habitat Stability Class`,`UCA Class`) %>% 
#   distinct() %>% 
#   group_by(`Habitat Stability Class`) %>% 
#   group_split() %>% 
#   map(~setNames(list(.x$`UCA Class`),.x$`Habitat Stability Class`[[1]])) %>% 
#   unlist(recursive = F)
# 
# site_grp <- site_grp[c("Low","High")]
# names(site_grp) <- c(c("Low Stability","High Stability"))
# 
# writexl::write_xlsx(shade_tbl,"data/changetbl3.xlsx")
# 
# shade_tbl <- shade_tbl %>% 
#   select(-`Habitat Stability Class`)
# 
# shade_tbl_sym <- shade_tbl %>% 
#   mutate(across(everything(),
#                 ~case_when(
#                   .x == "Increasing Significant" ~ "++",
#                   .x == "Increasing Trend" ~ "+",
#                   .x == "No Change" ~ "0",
#                   .x == "Decreasing Trend" ~ "-",
#                   .x == "Decreasing Significant" ~ "--",
#                   is.na(.x) ~ "0",
#                   T ~ .x
#                 )))
# 
# gt_tbl3 <- gt::gt(shade_tbl_sym,rowname_col = "UCA Class") %>% 
#   gt::tab_stubhead(label = "UCA Class / Habitat Stability Class") %>% 
#   gt::tab_source_note(
#     source_note = gt::md(col_list_leg)
#   ) %>% 
#   gt::tab_style(
#     style = gt::cell_text(align = "center", v_align = "middle"),
#     locations = gt::cells_body(
#       columns = everything(),
#       rows = everything()
#     )
#   ) %>% 
#   gt::tab_style(
#     style = gt::cell_text(align = "center",
#                           v_align = "middle",
#                           stretch = "extra-condensed",
#                           size = "small"),
#     locations = gt::cells_column_labels()
#   ) %>% 
#   gt::tab_style(
#     style = gt::cell_text(align = "center",
#                           v_align = "middle",
#                           stretch = "extra-condensed",
#                           size = "small",
#                           transform="capitalize"),
#     locations = gt::cells_row_groups()
#   )
# 
# for (i in names(col_list)){
#   for (ii in 1:ncol(shade_tbl)){
#     gt_tbl3 <- gt_tbl3 %>% 
#       gt::tab_style(
#         style = gt::cell_fill(color = col_list[[i]]),
#         locations = gt::cells_body(
#           columns = ii,
#           rows = which(shade_tbl[[ii]] == i)
#         )
#       )
#   }
# }
# 
# for (i in names(site_grp)){
#   gt_tbl3 <- gt_tbl3 %>% gt::tab_row_group(label = i,rows = site_grp[[i]])
# }
# 
# for (i in names(ep_grp)){
#   gt_tbl3 <- gt_tbl3 %>% gt::tab_spanner(label = i,columns  = as.character(ep_grp[[i]]))
# }

gt_tbl4 <- gt::gt_group(
  gt_tbl1,
  gt_tbl2,
  #gt_tbl3
)

saveRDS(gt_tbl4,"data/summarytbl2.rds")

6 Biological Correlations

source("00_source.R")

mod_data_fin <- readRDS("data/results_out.rds")

fish_ep <- read_csv("data/fish_endpoints.csv") %>% 
  select(SiteCode,Year,Endpoint,Value) %>% 
  pivot_wider(names_from = Endpoint, values_from = Value) %>% 
  select(SiteCode,Year, any_of(sel_fish_ep))

cor_list <- list()
fish_plot_list <- list()

# TRCA Data ---------------------------------------------------------------

cmb_data_TRCA <- mod_data_fin %>% 
  filter(TRCA) %>% 
  select(Endpoint, data) %>% 
  unnest(data) %>% 
  pivot_wider(names_from = Endpoint,values_from = value) %>% 
  left_join(fish_ep,by = join_by(SiteCode, Year))

hab_data_TRCA <- cmb_data_TRCA %>% 
  select(any_of(unlist(sapply(met_cols,names),use.names=F))#,
         # any_of(c("Upstream Catchment Area","Reach Slope",
         #          "Upstream catchment stress index","Upstream catchment developed (%)","Upstream catchment agricultural (%)"))
  )

fish_data_TRCA <- cmb_data_TRCA %>% 
  select(any_of(names(sel_fish_ep)))

cor_list$TRCA <- heatmap_fn(hab_data_TRCA,
                            fish_data_TRCA,
                            alpha=0.1,
                            sig_est=0.6,
                            method="spearman",
                            scale_grad=ggplot2::scale_fill_gradient2(low = scales::muted(rgb(red=0/255,green=102/255,blue=255/255)),
                                                                     mid = "white",
                                                                     high = scales::muted(rgb(red=0/255,green=176/255,blue=80/255)),
                                                                     limits=c(-1,1))) 

fish_plot_list$TRCA <- cor(fish_data_TRCA,hab_data_TRCA,method="spearman",use="pairwise.complete.obs") %>% 
  as.data.frame() %>% 
  rownames_to_column("var1") %>% 
  pivot_longer(c(everything(),-var1),names_to = "var2") %>% 
  mutate(abs_cor = abs(value)) %>% 
  filter(abs_cor>0.3)%>% 
  arrange(desc(abs_cor)) %>% 
  mutate(plot = map2(
    var1,
    var2,
    ~tibble(
      SiteCode = cmb_data_TRCA$StreamCode,
      SampleDate = cmb_data_TRCA$SampleDate,
      var2 = cmb_data_TRCA[[.y]],
      var1 = fish_data_TRCA[[.x]]
    ) %>% 
      ggplot(aes(x=var2,y=var1))+
      geom_point()+
      geom_smooth(se=F)+
      ylab(.x)+
      xlab(.y)+
      theme_bw()
    
  )) %>% 
  mutate(plot = map2(plot,value,
                     ~.x + labs(
                       #title = "TRCA Data Only",
                       subtitle = paste0("Spearman Correlation = ",scales::label_number(0.001)(.y)))
  ))

# All Data ----------------------------------------------------------------

cmb_data_all <- mod_data_fin %>% 
  filter(!TRCA) %>% 
  select(Endpoint, data) %>% 
  unnest(data) %>% 
  pivot_wider(names_from = Endpoint,values_from = value) %>% 
  left_join(fish_ep,by = join_by(SiteCode, Year))

hab_data_all <- cmb_data_all %>% 
  select(any_of(unlist(sapply(met_cols,names),use.names=F))#,
         # any_of(c("Upstream Catchment Area","Reach Slope",
         #          "Upstream catchment stress index","Upstream catchment developed (%)","Upstream catchment agricultural (%)"))
  )

fish_data_all <- cmb_data_all %>% 
  select(any_of(names(sel_fish_ep)))

cor_list$All <- heatmap_fn(hab_data_all,
                            fish_data_all,
                            alpha=0.1,
                            sig_est=0.6,
                            method="spearman",
                            scale_grad=ggplot2::scale_fill_gradient2(low = scales::muted(rgb(red=0/255,green=102/255,blue=255/255)),
                                                                     mid = "white",
                                                                     high = scales::muted(rgb(red=0/255,green=176/255,blue=80/255)),
                                                                     limits=c(-1,1))) 


fish_plot_list$All <- cor(fish_data_all,hab_data_all,method="spearman",use="pairwise.complete.obs") %>% 
  as.data.frame() %>% 
  rownames_to_column("var1") %>% 
  pivot_longer(c(everything(),-var1),names_to = "var2") %>% 
  mutate(abs_cor = abs(value)) %>% 
  filter(abs_cor>0.3)%>% 
  arrange(desc(abs_cor)) %>% 
  mutate(plot = map2(
    var1,
    var2,
    ~tibble(
      SiteCode = cmb_data_all$StreamCode,
      SampleDate = cmb_data_all$SampleDate,
      var2 = hab_data_all[[.y]],
      var1 = fish_data_all[[.x]]
    ) %>% 
      ggplot(aes(x=var2,y=var1))+
      geom_point()+
      geom_smooth(se=F)+
      ylab(.x)+
      xlab(.y)+
      theme_bw()
    
  )) %>% 
  mutate(plot = map2(plot,value,
                     ~.x + labs(
                       #title = "All Data",
                       subtitle = paste0("Spearman Correlation = ",scales::label_number(0.001)(.y)))
  ))

saveRDS(fish_plot_list,"data/biofish_plots.rds")
saveRDS(cor_list,"data/biocor_plots.rds")

Raw with Model Predictions Data

These figures show the raw data for each endpoint over time, grouped by watershed. The grey bars represent the 90 percent prediction interval for each observation.

Morphology

Width - Mean

Depth - Mean

Width:Depth Ratio

Usable Area

Substrate Sizes

D50

D85

D50Max

Indices

Sediment Sorting Index

Bank Stability - Mean

Habitat Homogeneity

Percent

Filamentous Algae Percent

Non-Filamentous Algae Percent

Moss Percent

Cover

Wood Cover - Unembedded

Wood Cover - Total

Round Rock Cover - Unembedded

Round Rock Cover - Total

Flat Rock Cover - Unembedded

Flat Rock Cover - Total

Results - Model Summaries

Model Summaries

These tables show the model coefficients, estimated standard error, and significance for fixed effects terms. Random effects are summarized as variance and correlations within groups. Variance measures how much the random effects vary across different groups or levels in the hierarchy.Random effects models assume that observations within the same group are correlated. This correlation arises because these observations share common characteristics. Summary measures are provided below the model results.

All Data

effect			Width - Mean		Depth - Mean		Width:Depth Ratio		Usable Area		D50		D85		D50Max		Sediment Sorting Index		Bank Stability - Mean		Habitat Homogeneity		Filamentous Algae Percent		Non-Filamentous Algae Percent		Moss Percent		Wood Cover - Unembedded		Wood Cover - Total		Round Rock Cover - Unembedded		Round Rock Cover - Total		Flat Rock Cover - Unembedded		Flat Rock Cover - Total
			Est.	S.E.	Est.	S.E.	Est.	S.E.	Est.	S.E.	Est.	S.E.	Est.	S.E.	Est.	S.E.	Est.	S.E.	Est.	S.E.	Est.	S.E.	Est.	S.E.	Est.	S.E.	Est.	S.E.	Est.	S.E.	Est.	S.E.	Est.	S.E.	Est.	S.E.	Est.	S.E.	Est.	S.E.
fixed cond	(Intercept)		0.943***	0.163	5.390***	0.143	2.482***	0.186	4.731***	0.203	2.761***	0.269	3.652***	0.287	3.271***	0.308	5.140***	0.188	0.120	0.142	-2.491***	0.328	-1.923**	0.649	-2.696***	0.718	-2.275***	0.523	-2.537***	0.283	-2.430***	0.304	-2.550***	0.391	-2.512***	0.425	-3.328***	0.452	-3.081***	0.479
	log1p(`Reach Slope`)		0.005	0.187	-0.593***	0.119	0.540**	0.169	-0.190	0.184	0.304+	0.168	0.932***	0.219	1.133***	0.250	-0.069	0.136	0.219*	0.106	0.441*	0.224	0.646+	0.356	0.942**	0.337	0.675*	0.330	-0.011	0.226	0.046	0.231	0.835**	0.281	1.220***	0.341	1.201***	0.354	1.574***	0.396
	Year_n		0.001	0.005	0.002	0.006	0.003	0.006	0.012+	0.006	-0.004	0.011	-0.006	0.010	-0.002	0.009	-0.002	0.008	0.004	0.006	0.029*	0.014	-0.011	0.034	0.001	0.032	0.003	0.022	0.022*	0.011	0.017	0.012	-0.010	0.016	-0.014	0.015	0.030+	0.017	0.015	0.017
	`UCA Class`Intermediate		0.707***	0.142	0.276+	0.147	0.503*	0.197	0.931***	0.220	-0.276	0.338	0.211	0.317	0.607+	0.313	-0.414+	0.222	-0.179	0.148	1.682***	0.361	0.282	0.703	0.484	0.758	0.876	0.558	0.081	0.306	0.149	0.333	0.865*	0.406	1.005*	0.444	0.926*	0.456	0.837+	0.472
	`UCA Class`Large		1.333***	0.264	0.335+	0.190	1.079***	0.281	1.719***	0.322	-0.158	0.366	0.602	0.400	0.718+	0.432	-0.634*	0.277	0.160	0.199	2.115***	0.458	0.279	0.783	0.476	0.885	0.903	0.833	-0.557	0.478	-0.545	0.520	0.846	0.545	1.158+	0.621	0.602	0.668	0.577	0.693
	`Stress Class`Medium		0.190	0.171	0.091	0.144	0.280	0.181	0.346+	0.200	0.181	0.291	-0.002	0.288	0.271	0.285	-0.249	0.209	-0.084	0.138	1.095**	0.352	-0.102	0.674	-0.140	0.721	1.566**	0.540	-0.457	0.332	-0.564	0.353	0.477	0.375	0.528	0.400	1.036*	0.411	0.915*	0.425
	`Stress Class`High		0.401**	0.128	0.154	0.132	0.323+	0.170	0.551**	0.190	-0.056	0.267	0.234	0.274	0.545*	0.276	-0.263	0.193	-0.065	0.129	0.653*	0.333	-0.320	0.661	0.216	0.685	1.188*	0.525	0.004	0.332	-0.089	0.354	0.555	0.364	0.733+	0.385	0.826*	0.395	0.769+	0.406
	Year_n × `UCA Class`Intermediate		0.001	0.006	0.000	0.007	0.003	0.009	-0.001	0.008	0.028+	0.017	0.010	0.014	0.003	0.012	0.010	0.011	0.008	0.008	-0.035*	0.017	0.018	0.037	0.040	0.037	-0.021	0.027	-0.005	0.015	-0.008	0.016	0.005	0.020	-0.003	0.019	-0.016	0.021	-0.003	0.021
	Year_n × `UCA Class`Large		-0.001	0.007	0.003	0.009	-0.012	0.013	-0.011	0.013	0.026	0.021	-0.002	0.018	-0.002	0.017	0.011	0.015	-0.013	0.011	-0.036	0.023	0.027	0.044	0.026	0.048	-0.042	0.054	0.002	0.030	0.003	0.032	-0.008	0.031	-0.016	0.029	-0.007	0.033	0.005	0.032
	Year_n × `Stress Class`Medium		0.002	0.006	-0.003	0.007	-0.004	0.008	-0.007	0.008	0.006	0.014	0.020	0.013	0.008	0.011	-0.000	0.010	0.006	0.007	-0.032+	0.016	0.054	0.035	0.063+	0.034	-0.054*	0.025	-0.005	0.015	0.001	0.016	0.013	0.018	0.009	0.017	-0.015	0.018	-0.004	0.018
	Year_n × `Stress Class`High		0.004	0.005	0.003	0.006	-0.004	0.007	-0.001	0.007	0.020	0.013	0.016	0.012	0.006	0.010	-0.006	0.010	-0.003	0.007	-0.014	0.016	0.065+	0.035	0.052	0.033	-0.039	0.026	-0.024	0.016	-0.022	0.017	0.022	0.017	0.010	0.016	0.006	0.017	0.012	0.017
	`UCA Class`Intermediate × `Stress Class`Medium		-0.334+	0.196	-0.183	0.184	-0.317	0.246	-0.499+	0.277	0.313	0.403	0.517	0.380	-0.137	0.375	0.235	0.277	0.140	0.185	-2.186***	0.451	-0.160	0.793	0.117	0.860	-1.009	0.676	0.082	0.481	0.002	0.509	-0.738	0.480	-0.792	0.533	-0.601	0.541	-0.304	0.564
	`UCA Class`Large × `Stress Class`Medium		-0.598*	0.269	-0.248	0.217	-0.122	0.314	-0.480	0.358	0.256	0.420	0.187	0.445	0.131	0.467	0.556+	0.321	-0.041	0.226	-2.030***	0.526	0.208	0.863	0.199	0.963	0.064	0.919	0.728	0.641	0.848	0.693	-0.696	0.599	-0.871	0.679	0.248	0.718	0.577	0.747
	`UCA Class`Intermediate × `Stress Class`High		-0.136	0.097	-0.151	0.188	-0.083	0.264	-0.244	0.300	0.285	0.438	-0.334	0.404	-0.692+	0.402	0.313	0.290	-0.050	0.193	-1.716***	0.476	-0.235	0.807	-0.337	0.845	-0.358	0.709	0.016	0.672	-0.289	0.669	-0.921+	0.494	-0.891	0.557	-0.887	0.567	-0.682	0.594
	`UCA Class`Large × `Stress Class`High		-0.083	0.276	-0.024	0.226	-0.229	0.338	-0.258	0.388	0.430	0.434	0.075	0.478	0.116	0.502	0.390	0.344	-0.066	0.242	-1.839**	0.562	-0.025	0.892	-0.056	0.963	-0.007	0.960	-0.066	1.388	-0.232	1.181	-1.597*	0.645	-1.945**	0.734	0.315	0.757	0.388	0.790
	Year_n × `UCA Class`Intermediate × `Stress Class`Medium		0.001	0.007	0.007	0.009	-0.007	0.011	0.002	0.011	-0.021	0.021	-0.032+	0.017	-0.011	0.014	-0.005	0.014	-0.007	0.010	0.058**	0.022	-0.022	0.041	-0.058	0.042	0.031	0.033	0.000	0.023	0.006	0.024	-0.009	0.023	-0.001	0.022	0.011	0.024	-0.007	0.024
	Year_n × `UCA Class`Large × `Stress Class`Medium		0.006	0.008	0.006	0.010	0.012	0.015	0.019	0.014	-0.019	0.024	-0.016	0.020	-0.004	0.019	-0.019	0.017	0.009	0.013	0.049+	0.026	-0.019	0.047	-0.020	0.052	0.024	0.058	-0.009	0.037	-0.017	0.039	0.010	0.033	0.019	0.031	0.012	0.036	-0.001	0.034
	Year_n × `UCA Class`Intermediate × `Stress Class`High		-0.001	0.007	0.003	0.009	-0.005	0.012	0.002	0.011	-0.016	0.023	0.018	0.017	0.022	0.015	-0.006	0.015	0.001	0.010	0.050*	0.023	0.005	0.042	-0.020	0.042	0.019	0.035	-0.003	0.032	0.010	0.032	0.017	0.023	0.026	0.022	0.032	0.025	0.020	0.025
	Year_n × `UCA Class`Large × `Stress Class`High		0.002	0.008	-0.006	0.011	0.018	0.016	0.002	0.015	-0.019	0.025	-0.011	0.021	-0.016	0.020	-0.013	0.018	0.012	0.014	0.057*	0.028	-0.027	0.048	-0.018	0.052	0.014	0.060	0.010	0.068	0.016	0.060	0.029	0.035	0.046	0.033	-0.008	0.037	-0.007	0.036
ran_pars cond	sd__(Intercept)	`UCA Class`:`Stress Class`:SiteCode	0.567		0.267		0.453		0.549		0.305		0.511		0.653		0.262		0.273		0.526		0.491		0.224		0.746		0.256		0.412		0.519		0.781		0.753		0.865
		Year_n	0.035		0.071		0.060		0.029		0.059		0.042		0.050		0.065		0.101		0.098		0.270		0.212		0.228		0.075		0.095		0.180		0.193		0.235		0.266
	sd__Year_n	`UCA Class`:`Stress Class`:SiteCode	0.009		0.009		0.017		0.017		0.013		0.012		0.019		0.006		0.014		0.016		0.012		0.018		0.009		0.006		0.015		0.009		0.013		0.022		0.024
	cor__(Intercept).Year_n		-0.485		-0.246		-0.485		-0.603		-0.463		-0.439		-0.548		-0.213		-0.553		-0.492		1.000		1.000		-1.000		0.806		-0.477		0.260		-0.297		-0.155		-0.161
fixed zi	(Intercept)										4.663***	1.378	-0.795	3.450	-1.760	3.869	-6.867+	4.049			-0.409	1.281	2.160**	0.724	4.340***	0.759	-0.687	0.828	0.318	0.818	-0.079	0.830	3.788**	1.335	1.803	1.876	4.053*	1.815	0.683	4.285
	log1p(`Reach Slope`)										-6.763***	1.680	-7.184	4.961	-4.585	3.745	-3.338	4.338			-1.558	1.292	-0.318	0.654	-1.384*	0.615	-3.356***	0.856	0.997	0.799	1.518+	0.812	-7.046***	1.853	-8.234***	2.112	-7.026**	2.604	-8.490+	4.847
	Year_n										0.002	0.031	0.018	0.138	-0.135	0.137	0.032	0.186			-0.177**	0.061	-0.096***	0.028	-0.129***	0.030	0.159***	0.025	-0.117***	0.022	-0.119***	0.022	-0.017	0.036	0.059	0.055	-0.095	0.080	-0.124	0.140
	`UCA Class`Intermediate										-1.961**	0.759	-2.491	1.568	-1.998	1.662	-0.833	2.168			-1.286*	0.652	-1.033**	0.332	-1.290***	0.310	-1.753***	0.431	0.643	0.413	0.639	0.418	-2.310**	0.802	-2.262*	0.988	-2.936**	1.054	-2.623	1.728
	`UCA Class`Large										-4.961***	1.011	-4.961+	2.858	-5.159	3.277	-18.505	8556.698			-1.970*	0.771	-1.594***	0.418	-2.061***	0.389	-0.633	0.514	2.170***	0.524	2.412***	0.532	-2.458*	1.047	-2.718*	1.154	-3.132*	1.423	-3.713	2.490
	`Stress Class`Medium										-1.355+	0.737	-2.516+	1.409	-2.202	1.624	-2.845	2.474			0.794	0.654	-1.022**	0.342	-1.005**	0.309	-0.302	0.418	1.623***	0.422	1.684***	0.428	-2.408**	0.787	-2.973**	0.952	-3.526***	1.035	-3.065+	1.765
	`Stress Class`High										-2.370**	0.847	-3.279+	1.744	-2.251	1.812	-2.522	2.391			0.156	0.702	-1.291***	0.386	-2.241***	0.372	-0.277	0.462	2.488***	0.471	2.579***	0.480	-2.771**	0.905	-3.396**	1.074	-3.796***	1.129	-3.657+	1.976
ran_pars zi	sd__(Intercept)	`UCA Class`:`Stress Class`:SiteCode									2.755		6.017		5.897		12.833				0.762		0.898		1.286		2.491		2.291		2.572		2.815		5.200		3.530		10.250
		Year_n									0.710		0.796		0.678		0.001				0.323		0.821		0.865		0.591		0.426		0.420		0.195		0.391		0.777		1.398
	sd__Year_n	`UCA Class`:`Stress Class`:SiteCode									0.074		0.010		0.430		0.660				0.098		0.029		0.065		0.082		0.055		0.055		0.000		0.118		0.160		0.426
	cor__(Intercept).Year_n										-0.038		0.997		0.106		-0.840				1.000		0.557		-0.608		-0.767		-0.795		-0.931		0.997		-0.819		-0.104		-0.602
	Num.Obs.		1168		1168		1168		1157		1059		1078		1091		1166		1177		1168		1168		1168		1168		1168		1168		1168		1168		1168		1168
	RMSE		5.01		54.70		10.91		82.94		12.47		72.65		37.87		52.13		0.08		0.11		0.16		0.16		0.14		0.04		0.04		0.09		0.09		0.09		0.09
	Pseudo R2		0.771		0.822		0.851		0.937		0.367		0.728		0.783		0.331		0.638		0.547		0.485		0.386		0.317		0.180		0.203		0.602		0.707		0.763		0.809
+ p < 0.1, * p < 0.05, ** p < 0.01, *** p < 0.001
zi = zero-inflation; cond = conditional; fixed = fixed effects parameter; ran_pars = random effects parameters;

Model Residuals

Results of the DHARMa package which uses a simulation-based approach to create readily interpretable scaled (quantile) residuals for fitted generalized linear (mixed) models. Left panel is a qq-plot to detect overall deviations from the expected distribution, with added tests for correct distribution (KS test), dispersion and outliers. Note that outliers are values that are defined as values outside the simulation envelope, not in terms of a particular quantile. Right panel is a plot of the residuals against the predicted value. Simulation outliers (data points that are outside the range of simulated values) are highlighted as red stars.

Morphology

Width - Mean

Depth - Mean

Width:Depth Ratio

Usable Area

Substrate Sizes

D50

D85

D50Max

Indices

Sediment Sorting Index

Bank Stability - Mean

Habitat Homogeneity

Percent

Filamentous Algae Percent

Non-Filamentous Algae Percent

Moss Percent

Cover

Wood Cover - Unembedded

Wood Cover - Total

Round Rock Cover - Unembedded

Round Rock Cover - Total

Flat Rock Cover - Unembedded

Flat Rock Cover - Total

Model Fit

These figures show 1:1 plots of observed and model fitted values. The plots are shown on transformed scales. Endpoints were model fitted values fall more closely to the 1:1 line better fit the raw data. A good fitting model means the endpoint align with the distributional family in both mean and variance, that means of monitoring locations are normally distributed (on the transformed scale), and the overdispersion and zero-inflation terms are able to describe the attributes of the raw data. A pseudo R2 value is shown as represented by the squared spearman correlation coefficient of observed vs predicted values.

Morphology

Width - Mean

Depth - Mean

Width:Depth Ratio

Usable Area

Substrate Sizes

D50

D85

D50Max

Indices

Sediment Sorting Index

Bank Stability - Mean

Habitat Homogeneity

Percent

Filamentous Algae Percent

Non-Filamentous Algae Percent

Moss Percent

Cover

Wood Cover - Unembedded

Wood Cover - Total

Round Rock Cover - Unembedded

Round Rock Cover - Total

Flat Rock Cover - Unembedded

Flat Rock Cover - Total

Results - Spatial Figures

Overall

This figure shows the global means and associated error bars. The Overall mean reflects the patterns across all monitoring locations, whereas the UCA Class and Stress Class, each reflect the patterns across locations assigned to those levels. For more detailed breakdowns across combinations of these 3 location classes, see the UCA and Stress tab.

UCA and Stress

Endpoints by Monitoring Locations

Morphology

Width - Mean

Depth - Mean

Width:Depth Ratio

Usable Area

Substrate Sizes

D50

D85

D50Max

Indices

Sediment Sorting Index

Bank Stability - Mean

Habitat Homogeneity

Percent

Filamentous Algae Percent

Non-Filamentous Algae Percent

Moss Percent

Cover

Wood Cover - Unembedded

Wood Cover - Total

Round Rock Cover - Unembedded

Round Rock Cover - Total

Flat Rock Cover - Unembedded

Flat Rock Cover - Total

Spatial Correlations

Heatmap of spearman correlation coefficients among habitat endpoints (including Aquatic Ecosystem Classification variables), with green and blue indicating positive and negative relationships respectively. Strong Statistically significant correlations (rho > 0.6 & p < 0.1) are shown with red outlines. Moderate statistically significant correlations (rho > 0.4 & p < 0.1) are shown with orange outlines. Endpoints are arranged according to a Ward clustering.

Results - Temporal Figures

Overall

This figure shows the global temporal trend coefficients and associated error bars. The Overall trend reflects the patterns across all monitoring locations, whereas the UCA Class and Stress Class each reflect the patterns across locations assigned to those levels. For more detailed breakdowns across combinations of these 3 location classes, see the UCA and Stress tab.

UCA and Stress

Endpoints by Monitoring Locations

Morphology

Width - Mean

Depth - Mean

Width:Depth Ratio

Usable Area

Substrate Sizes

D50

D85

D50Max

Indices

Sediment Sorting Index

Bank Stability - Mean

Habitat Homogeneity

Percent

Filamentous Algae Percent

Non-Filamentous Algae Percent

Moss Percent

Cover

Wood Cover - Unembedded

Wood Cover - Total

Round Rock Cover - Unembedded

Round Rock Cover - Total

Flat Rock Cover - Unembedded

Flat Rock Cover - Total

Trend Correlation

Heatmap and cluster diagram of spearman correlations among temporal trends in habitat endpoints. Trends were coded as +2, +1, 0, -1, and -2 corresponding with ++, +, 0, -, and -- from the Trends Table. Strong positive correlations among habitat endpoints indicate shared increasing or decreasing temporal trends across multiple sites, while strong negative correlations indicate opposing increasing and decreasing temporal trends across multiple sites. Strong Statistically significant correlations (rho > 0.6 & p < 0.1) are shown with red outlines. Moderate statistically significant correlations (rho > 0.4 & p < 0.1) are shown with orange outlines. Endpoints are arranged according to a Ward clustering.

Results - Temporal Tables

Trends by Stream Order and Stress

This table shows concurrent patterns of temporal trends across all habitat endpoints by stream order and stress class. The ++ and -- indicate that the 90 percent confidence interval of the slope coefficient does not intersect 0, suggesting a significant increasing and decreasing patterns, respectively. The + and - indicate that the 50 percent confidence interval of the of the slope coefficient does not intersect 0, suggesting increasing and decreasing non-significant trends.

Overall Trends	Morphology				Substrate Sizes			Indices			Percent			Cover
	Width - Mean	Depth - Mean	Width:Depth Ratio	Usable Area	D50	D85	D50Max	Sediment Sorting Index	Bank Stability - Mean	Habitat Homogeneity	Filamentous Algae Percent	Non-Filamentous Algae Percent	Moss Percent	Wood Cover - Unembedded	Wood Cover - Total	Round Rock Cover - Unembedded	Round Rock Cover - Total	Flat Rock Cover - Unembedded	Flat Rock Cover - Total
Global Trend	++	+	0	++	++	+	0	0	+	++	++	++	--	+	+	+	0	++	++
Stress Class
Low	0	+	0	+	+	0	0	+	0	0	0	+	-	+	+	-	-	+	+
Medium	++	+	0	++	+	0	0	0	++	+	++	++	--	+	+	0	0	+	+
High	++	+	0	++	++	++	+	-	+	++	++	++	--	0	0	++	+	++	++
UCA Class
Small	+	0	0	++	+	+	0	-	+	++	+	++	--	+	+	0	-	++	+
Intermediate	+	+	0	++	++	++	++	+	++	++	++	++	--	0	+	+	0	++	+
Large	+	+	0	+	++	-	-	0	0	+	++	++	--	0	0	0	0	+	+
+ + = Increasing Significant | + = Increasing Trend | 0 = No Change | - = Decreasing Trend | - - = Decreasing Significant

UCA Class / Stress Class	Morphology				Substrate Sizes			Indices			Percent			Cover
	Width - Mean	Depth - Mean	Width:Depth Ratio	Usable Area	D50	D85	D50Max	Sediment Sorting Index	Bank Stability - Mean	Habitat Homogeneity	Filamentous Algae Percent	Non-Filamentous Algae Percent	Moss Percent	Wood Cover - Unembedded	Wood Cover - Total	Round Rock Cover - Unembedded	Round Rock Cover - Total	Flat Rock Cover - Unembedded	Flat Rock Cover - Total
High Stress
UCA Small	+	+	0	++	++	+	+	-	0	+	++	++	--	0	0	+	0	++	++
UCA Intermediate	+	+	0	+	++	++	++	0	+	++	++	++	--	0	0	++	+	++	++
UCA Large	+	0	0	0	++	0	-	-	0	++	++	++	--	0	0	++	+	+	+
Medium Stress
UCA Small	+	0	0	+	0	+	+	0	+	0	++	++	--	+	+	0	0	+	+
UCA Intermediate	+	+	-	+	+	-	0	0	++	++	++	++	--	+	+	0	-	+	0
UCA Large	++	+	0	++	+	0	0	-	+	+	++	++	--	0	0	0	0	+	+
Low Stress
UCA Small	0	0	0	+	0	0	0	0	0	++	0	0	0	++	+	0	-	+	+
UCA Intermediate	0	0	+	+	+	0	0	+	+	0	0	+	-	+	+	0	-	+	+
UCA Large	0	0	-	0	+	0	0	+	-	0	0	+	-	+	+	-	-	+	+
+ + = Increasing Significant | + = Increasing Trend | 0 = No Change | - = Decreasing Trend | - - = Decreasing Significant

Trends by Location

This table shows concurrent patterns of temporal trends across all habitat endpoints by monitoring locations. The ++ and -- indicate that the 90 percent confidence interval of the slope coefficient does not intersect 0, suggesting a significant increasing and decreasing patterns, respectively. The + and - indicate that the 50 percent confidence interval of the of the slope coefficient does not intersect 0, suggesting increasing and decreasing non-significant trends.

HUMBER
HUMBER	UCA Class	Morphology				Substrate Sizes			Indices			Percent			Cover
		Width - Mean	Depth - Mean	Width:Depth Ratio	Usable Area	D50	D85	D50Max	Sediment Sorting Index	Bank Stability - Mean	Habitat Homogeneity	Filamentous Algae Percent	Non-Filamentous Algae Percent	Moss Percent	Wood Cover - Unembedded	Wood Cover - Total	Round Rock Cover - Unembedded	Round Rock Cover - Total	Flat Rock Cover - Unembedded	Flat Rock Cover - Total
High Stess Class
HU008WM	Small	+	+	0	+	+	+	++	-	-	0	++	++	--	0	0	+	0	++	++
HU005WM	Small	0	+	-	+	+	+	0	0	-	+	++	++	--	0	0	0	0	+	+
HU006WM	Intermediate	-	0	0	-	+	++	++	0	0	++	++	++	-	0	0	++	+	++	++
Medium Stess Class
HU039WM	Small	0	0	-	0	0	+	0	0	+	0	++	++	--	+	+	0	0	0	0
HU016WM	Small	0	0	+	++	0	+	0	0	+	0	++	++	--	+	+	0	0	+	0
HU001WM	Intermediate	0	+	0	-	0	0	0	0	0	+	++	++	--	+	0	0	0	0	0
HU028WM	Intermediate	0	+	-	0	+	0	0	0	+	+	++	++	--	+	+	0	0	0	0
HU026WM	Intermediate	+	+	0	+	0	0	0	0	0	+	++	+	--	+	+	0	-	0	0
HU027WM	Intermediate	0	0	+	+	0	-	0	0	++	+	++	++	--	+	0	0	0	+	0
HU011WM	Intermediate	+	0	+	0	+	0	0	0	0	+	++	++	--	+	0	0	0	+	+
HU015WM	Intermediate	0	0	0	0	0	-	0	0	+	+	++	++	--	+	+	0	0	+	+
HU014WM	Intermediate	0	0	0	0	0	-	0	0	0	+	++	++	--	+	+	0	0	0	0
HU022WM	Large	0	+	0	+	0	0	0	-	+	+	++	++	--	0	0	0	0	0	0
HU024WM	Large	0	+	0	++	0	0	0	0	0	0	+	++	--	+	0	0	0	0	0
HU010WM	Large	0	+	0	+	++	0	0	-	0	0	++	++	--	0	0	+	0	+	+
HU003WM	Large	++	0	0	+	0	-	0	-	-	+	++	++	--	0	0	0	0	++	+
HU007WM	Large	++	+	0	0	0	0	-	-	-	0	++	++	--	0	0	0	0	+	++
HU012WM	Large	0	+	0	-	0	0	0	0	-	0	++	++	--	0	0	0	0	+	+
HU013WM	Large	0	+	0	+	+	0	0	-	+	0	++	++	--	0	0	0	0	+	+
HU009WM	Large	+	0	0	+	0	0	0	-	+	+	++	++	--	0	0	0	0	0	0
Low Stess Class
HU032WM	Small	0	+	0	+	0	0	0	0	++	+	0	0	0	+	0	0	-	+	0
HU030WM	Small	0	0	+	+	0	0	0	0	-	+	0	0	0	++	+	-	-	+	+
HU035WM	Small	0	0	0	+	0	0	0	0	0	++	0	0	0	++	+	-	0	+	+
HU036WM	Intermediate	0	0	0	0	+	0	0	+	+	0	0	+	0	+	0	0	-	0	0
HU033WM	Intermediate	0	0	0	0	+	0	0	+	0	0	0	+	-	+	0	0	-	0	0
HU031WM	Intermediate	0	0	++	++	+	0	0	+	++	0	0	+	-	+	0	0	-	0	0
HU037WM	Intermediate	0	+	0	+	+	0	0	+	++	0	0	+	-	++	+	0	-	0	0
HU034WM	Intermediate	0	0	0	+	+	0	0	+	0	0	0	+	0	++	0	0	-	0	0
HU038WM	Intermediate	0	0	0	+	+	0	--	+	0	0	0	+	-	+	0	0	-	0	0
HU029WM	Large	0	+	--	0	0	0	-	0	-	0	0	+	-	+	+	-	-	0	0
HU002WM	Large	0	+	0	-	+	0	0	+	-	0	0	0	-	+	0	-	-	0	0
+ + = Increasing Significant | + = Increasing Trend | 0 = No Change | - = Decreasing Trend | - - = Decreasing Significant

CARRUTHERS
CARRUTHERS	UCA Class	Morphology				Substrate Sizes			Indices			Percent			Cover
		Width - Mean	Depth - Mean	Width:Depth Ratio	Usable Area	D50	D85	D50Max	Sediment Sorting Index	Bank Stability - Mean	Habitat Homogeneity	Filamentous Algae Percent	Non-Filamentous Algae Percent	Moss Percent	Wood Cover - Unembedded	Wood Cover - Total	Round Rock Cover - Unembedded	Round Rock Cover - Total	Flat Rock Cover - Unembedded	Flat Rock Cover - Total
Medium Stess Class
CC003WM	Small	0	0	0	+	0	+	-	0	0	0	++	++	--	+	+	0	0	0	0
CC002WM	Intermediate	++	++	-	++	0	0	-	+	+	+	++	++	--	+	+	0	0	0	0
CC001WM	Intermediate	+	0	0	++	0	0	0	+	0	+	++	+	--	+	+	0	0	0	0
+ + = Increasing Significant | + = Increasing Trend | 0 = No Change | - = Decreasing Trend | - - = Decreasing Significant

HIGHLAND
HIGHLAND	UCA Class	Morphology				Substrate Sizes			Indices			Percent			Cover
		Width - Mean	Depth - Mean	Width:Depth Ratio	Usable Area	D50	D85	D50Max	Sediment Sorting Index	Bank Stability - Mean	Habitat Homogeneity	Filamentous Algae Percent	Non-Filamentous Algae Percent	Moss Percent	Wood Cover - Unembedded	Wood Cover - Total	Round Rock Cover - Unembedded	Round Rock Cover - Total	Flat Rock Cover - Unembedded	Flat Rock Cover - Total
High Stess Class
HL001WM	Small	+	+	--	--	+	+	0	-	0	+	++	++	--	0	0	0	0	+	+
HL006WM	Small	0	0	+	0	+	0	0	-	+	+	++	++	--	0	0	+	0	++	++
HL008WM	Small	+	0	+	++	+	+	0	-	-	+	++	++	--	0	0	+	0	++	+
HL007WM	Small	0	0	--	0	+	0	0	-	+	+	++	++	--	0	0	0	0	+	0
HL009WM	Small	--	0	-	-	++	0	+	-	+	+	++	++	--	0	0	+	0	++	++
HL005WM	Intermediate	++	0	++	++	+	++	++	0	++	++	++	++	--	0	0	++	+	++	++
HL010WM	Intermediate	++	++	0	++	+	++	+	0	+	+	++	++	--	0	0	++	+	++	+
HL003WM	Intermediate	++	++	-	+	+	++	+	0	0	+	++	++	--	0	0	++	0	++	+
HL011WM	Intermediate	+	0	+	-	+	++	+	0	0	+	++	++	-	0	0	++	+	++	+
+ + = Increasing Significant | + = Increasing Trend | 0 = No Change | - = Decreasing Trend | - - = Decreasing Significant

DUFFINS
DUFFINS	UCA Class	Morphology				Substrate Sizes			Indices			Percent			Cover
		Width - Mean	Depth - Mean	Width:Depth Ratio	Usable Area	D50	D85	D50Max	Sediment Sorting Index	Bank Stability - Mean	Habitat Homogeneity	Filamentous Algae Percent	Non-Filamentous Algae Percent	Moss Percent	Wood Cover - Unembedded	Wood Cover - Total	Round Rock Cover - Unembedded	Round Rock Cover - Total	Flat Rock Cover - Unembedded	Flat Rock Cover - Total
High Stess Class
DF001WM	Small	--	+	0	0	+	+	0	-	-	0	++	+	--	0	0	+	0	+	+
DF002WM	Small	0	-	++	+	+	0	0	-	0	+	++	++	--	0	0	+	0	++	++
Medium Stess Class
DF018WM	Small	0	-	+	0	0	0	+	0	0	0	++	++	--	+	+	0	0	0	0
DF008WM	Small	+	0	+	++	0	+	0	0	0	0	+	++	--	++	+	0	0	0	0
DF005WM	Small	0	0	-	+	0	+	0	0	+	0	++	++	--	+	+	0	0	+	+
DF007WM	Small	+	0	+	+	0	+	0	0	0	0	++	++	--	+	+	0	0	0	0
DF014WM	Small	0	0	0	-	0	+	+	0	0	0	++	++	--	+	+	0	0	+	+
DF016WM	Small	0	0	0	0	0	+	0	0	+	0	++	++	--	+	+	0	0	+	+
DF020WMB	Intermediate	0	0	0	0	+	0	+	0	++	+	++	++	--	+	+	0	0	+	+
DF013WM	Intermediate	0	0	0	0	0	0	+	0	++	+	++	++	--	+	0	0	0	+	+
DF011WM	Large	0	0	0	+	0	0	0	-	++	+	++	++	--	0	0	0	0	+	+
DF010WM	Large	0	+	0	0	0	0	--	-	+	0	++	++	--	0	0	0	0	0	0
DF009WM	Large	++	+	0	++	0	0	-	0	0	0	++	++	--	0	0	0	0	+	0
DF003WM	Large	0	+	0	++	0	0	+	-	0	+	++	++	--	0	0	0	0	+	+
Low Stess Class
DF021WMB	Small	0	0	0	+	0	0	0	0	0	+	0	0	0	+	+	0	-	+	0
DF004WM	Small	0	0	+	+	0	0	0	0	0	+	0	0	0	+	0	0	0	+	+
09263_UP	Small	0	0	0	+	0	0	0	0	+	+	0	0	0	++	+	0	0	+	0
09263_DOWN	Small	0	0	0	0	0	0	0	0	+	+	0	0	0	++	+	0	0	+	+
DF019WM	Intermediate	0	0	+	0	+	0	0	+	+	0	+	++	-	+	0	0	-	0	0
DF017WM	Intermediate	-	0	0	0	+	0	+	+	0	0	0	+	0	+	+	-	-	0	0
DF015WM	Intermediate	0	0	+	+	+	0	0	0	+	0	0	++	-	+	0	0	-	+	+
DF012WM	Intermediate	0	0	0	0	++	+	+	0	+	0	+	++	-	+	0	0	0	+	+
DF006WM	Large	0	0	0	+	+	0	0	+	0	0	0	0	-	+	0	0	-	+	+
+ + = Increasing Significant | + = Increasing Trend | 0 = No Change | - = Decreasing Trend | - - = Decreasing Significant

ETOBICOKE
ETOBICOKE	UCA Class	Morphology				Substrate Sizes			Indices			Percent			Cover
		Width - Mean	Depth - Mean	Width:Depth Ratio	Usable Area	D50	D85	D50Max	Sediment Sorting Index	Bank Stability - Mean	Habitat Homogeneity	Filamentous Algae Percent	Non-Filamentous Algae Percent	Moss Percent	Wood Cover - Unembedded	Wood Cover - Total	Round Rock Cover - Unembedded	Round Rock Cover - Total	Flat Rock Cover - Unembedded	Flat Rock Cover - Total
High Stess Class
EC010WM	Small	+	+	+	++	+	+	0	-	0	+	++	++	--	0	0	+	0	++	+
EC006WM	Small	+	+	0	++	0	0	0	0	+	+	++	++	--	0	0	+	0	0	0
EC003WM	Small	+	+	0	++	0	0	--	0	-	+	++	++	--	0	0	0	0	0	-
EC009WM	Intermediate	-	+	-	-	+	++	++	0	+	+	++	++	--	0	0	++	+	++	++
EC005WM	Large	+	0	0	0	+	0	0	0	+	++	++	++	--	0	0	++	++	0	0
EC004WM	Large	++	-	++	0	++	0	--	0	-	++	++	++	--	0	0	+	+	+	+
EC002WM	Large	--	+	--	--	+	0	-	-	0	++	++	++	--	0	0	++	+	+	++
EC001WM	Large	0	++	-	0	+	0	0	-	0	++	++	++	--	0	0	+	+	++	+
Medium Stess Class
EC012WM	Intermediate	0	0	0	0	+	-	-	0	+	+	++	++	--	+	+	0	0	0	0
EC008WM	Large	+	0	0	+	0	-	0	-	+	+	++	++	--	0	0	0	0	+	+
EC007WM	Large	+	+	0	+	+	--	--	-	0	0	++	++	--	0	0	0	0	+	+
EC011WM	Large	0	0	0	0	+	0	0	-	0	+	++	++	--	0	0	0	0	+	+
+ + = Increasing Significant | + = Increasing Trend | 0 = No Change | - = Decreasing Trend | - - = Decreasing Significant

FRENCHMANS BAY
FRENCHMANS BAY	UCA Class	Morphology				Substrate Sizes			Indices			Percent			Cover
		Width - Mean	Depth - Mean	Width:Depth Ratio	Usable Area	D50	D85	D50Max	Sediment Sorting Index	Bank Stability - Mean	Habitat Homogeneity	Filamentous Algae Percent	Non-Filamentous Algae Percent	Moss Percent	Wood Cover - Unembedded	Wood Cover - Total	Round Rock Cover - Unembedded	Round Rock Cover - Total	Flat Rock Cover - Unembedded	Flat Rock Cover - Total
High Stess Class
FB001WM	Small	+	0	++	++	+	+	0	0	0	0	++	++	--	0	0	0	0	+	+
FB002WM	Small	+	0	-	+	+	+	+	-	0	0	++	++	-	0	0	+	0	+	+
FB004WM	Small	++	++	0	+	+	+	0	-	0	0	++	++	--	0	0	+	0	+	0
+ + = Increasing Significant | + = Increasing Trend | 0 = No Change | - = Decreasing Trend | - - = Decreasing Significant

DON
DON	UCA Class	Morphology				Substrate Sizes			Indices			Percent			Cover
		Width - Mean	Depth - Mean	Width:Depth Ratio	Usable Area	D50	D85	D50Max	Sediment Sorting Index	Bank Stability - Mean	Habitat Homogeneity	Filamentous Algae Percent	Non-Filamentous Algae Percent	Moss Percent	Wood Cover - Unembedded	Wood Cover - Total	Round Rock Cover - Unembedded	Round Rock Cover - Total	Flat Rock Cover - Unembedded	Flat Rock Cover - Total
High Stess Class
DN011WM	Small	+	+	0	+	0	0	--	-	0	+	++	++	--	0	0	0	-	++	+
DN010WMB	Small	+	+	0	+	+	0	-	-	0	+	++	++	--	0	0	+	0	+	+
DN009WMB	Small	0	0	0	0	+	+	+	-	0	+	++	++	--	0	0	+	0	++	++
DN006WMB	Small	0	++	-	0	0	0	0	-	+	+	++	++	--	0	0	+	0	++	+
DN007WM	Small	0	0	+	0	+	0	0	-	-	0	++	++	-	0	0	+	0	++	++
DN013WM	Small	+	0	+	++	+	+	++	-	-	+	++	++	--	0	0	+	0	++	++
DN014WM	Small	0	++	--	0	+	++	++	-	0	+	++	++	--	0	0	+	0	+	+
DN020WM	Small	0	0	0	0	+	+	0	-	0	+	++	++	--	0	0	0	0	+	+
DN018WM	Small	+	0	+	0	+	+	0	0	0	0	++	++	--	0	0	0	0	+	+
DN008WM	Intermediate	-	0	0	0	+	++	+	0	+	+	++	++	-	0	0	++	+	++	+
DN003WM	Intermediate	+	++	-	+	++	++	++	0	0	++	++	++	--	0	0	++	+	++	++
DN002WM	Intermediate	0	0	0	+	+	++	+	0	+	+	++	++	--	0	0	++	+	++	++
DN016WM	Intermediate	+	++	--	+	++	++	+	0	0	+	++	++	--	0	0	++	+	+	+
DN001WM	Large	++	0	+	+	+	0	0	-	0	+	++	+	--	0	0	+	+	0	0
DN005WM	Large	0	0	0	0	+	0	0	0	0	++	++	++	--	0	0	++	+	+	+
DN012WM	Large	+	0	+	+	+	0	-	-	0	+	++	++	--	0	0	+	+	0	0
DN015WM	Large	0	0	0	0	+	0	+	-	0	+	++	++	--	0	0	++	+	+	+
Medium Stess Class
DN023WM	Small	0	0	-	0	0	0	0	0	+	0	++	++	--	++	+	0	0	+	+
DN019WM	Small	0	0	+	+	0	+	+	0	+	0	++	++	--	+	+	0	0	0	0
DN021WM	Intermediate	+	+	0	+	+	-	0	0	+	+	++	++	--	+	+	0	0	0	0
Low Stess Class
DN022WM	Small	0	0	0	0	0	0	0	0	--	+	0	0	0	++	+	0	-	+	0
+ + = Increasing Significant | + = Increasing Trend | 0 = No Change | - = Decreasing Trend | - - = Decreasing Significant

ROUGE
ROUGE	UCA Class	Morphology				Substrate Sizes			Indices			Percent			Cover
		Width - Mean	Depth - Mean	Width:Depth Ratio	Usable Area	D50	D85	D50Max	Sediment Sorting Index	Bank Stability - Mean	Habitat Homogeneity	Filamentous Algae Percent	Non-Filamentous Algae Percent	Moss Percent	Wood Cover - Unembedded	Wood Cover - Total	Round Rock Cover - Unembedded	Round Rock Cover - Total	Flat Rock Cover - Unembedded	Flat Rock Cover - Total
High Stess Class
RG004WM	Small	0	0	0	0	+	0	0	-	+	+	++	++	--	0	0	+	0	+	0
RG005WM	Small	0	+	-	0	+	0	0	0	0	+	++	++	--	0	0	+	0	+	0
RG003WM	Small	+	+	0	++	0	0	0	-	0	+	++	++	--	0	0	+	0	+	0
RG017WM	Small	+	0	+	++	+	0	0	0	0	+	++	++	--	0	0	+	0	0	0
RG024WM	Small	+	0	+	+	+	+	+	-	-	+	++	++	--	0	0	+	0	+	+
Medium Stess Class
RG010WM	Small	0	0	0	0	+	+	0	0	+	0	++	++	--	+	+	+	0	+	+
RG022WM	Small	+	0	+	+	0	+	+	0	0	0	++	++	--	+	+	0	0	+	+
RG023WM	Small	0	0	+	0	0	+	0	0	+	0	++	++	--	+	+	0	0	+	+
RG011WM	Intermediate	+	+	0	+	0	-	-	0	0	+	+	++	--	+	0	0	0	0	0
RG018WM	Intermediate	0	0	0	+	+	0	0	0	+	+	++	++	--	+	+	0	-	+	+
RG015WM	Intermediate	+	0	0	+	+	-	0	0	+	+	++	++	--	+	+	0	0	0	0
RG016WM	Intermediate	0	+	0	0	0	-	-	0	+	+	++	++	--	+	0	0	-	+	+
RG019WM	Intermediate	+	++	--	+	0	0	0	0	+	+	+	+	--	+	+	0	0	0	0
RG020WM	Intermediate	+	++	--	+	0	-	0	0	0	+	++	++	--	+	+	0	0	0	0
RG021WM	Intermediate	0	0	-	0	0	0	0	0	0	+	++	++	--	+	+	0	-	0	0
RG007WM	Large	++	+	0	++	+	0	-	-	++	0	++	++	--	0	0	0	0	0	0
RG001WM	Large	+	+	0	++	0	0	++	0	++	+	++	++	--	0	0	0	0	0	0
RG006WM	Large	+	0	0	0	0	0	0	-	0	0	++	++	--	0	0	0	0	0	0
RG009WM	Large	-	++	-	++	+	0	++	-	0	+	++	++	--	0	0	0	0	+	+
RG002WM	Large	+	+	-	+	0	-	+	-	0	0	++	++	--	0	0	0	0	+	+
Low Stess Class
RG026WM	Small	-	0	++	++	0	0	0	0	0	+	0	0	0	++	+	-	-	+	0
RG027WM	Small	0	0	0	+	0	0	0	0	0	+	0	0	0	++	+	-	0	+	0
RG025WM	Small	0	0	-	0	0	0	-	0	0	+	0	0	0	+	++	0	-	+	0
+ + = Increasing Significant | + = Increasing Trend | 0 = No Change | - = Decreasing Trend | - - = Decreasing Significant

MIMICO
MIMICO	UCA Class	Morphology				Substrate Sizes			Indices			Percent			Cover
		Width - Mean	Depth - Mean	Width:Depth Ratio	Usable Area	D50	D85	D50Max	Sediment Sorting Index	Bank Stability - Mean	Habitat Homogeneity	Filamentous Algae Percent	Non-Filamentous Algae Percent	Moss Percent	Wood Cover - Unembedded	Wood Cover - Total	Round Rock Cover - Unembedded	Round Rock Cover - Total	Flat Rock Cover - Unembedded	Flat Rock Cover - Total
High Stess Class
MM005WM	Small	+	+	0	+	+	++	+	-	0	+	++	++	--	0	0	+	0	++	++
MM004WM	Intermediate	+	++	-	+	+	++	+	0	0	+	++	++	-	0	0	++	+	+	+
+ + = Increasing Significant | + = Increasing Trend | 0 = No Change | - = Decreasing Trend | - - = Decreasing Significant

PETTICOAT
PETTICOAT	UCA Class	Morphology				Substrate Sizes			Indices			Percent			Cover
		Width - Mean	Depth - Mean	Width:Depth Ratio	Usable Area	D50	D85	D50Max	Sediment Sorting Index	Bank Stability - Mean	Habitat Homogeneity	Filamentous Algae Percent	Non-Filamentous Algae Percent	Moss Percent	Wood Cover - Unembedded	Wood Cover - Total	Round Rock Cover - Unembedded	Round Rock Cover - Total	Flat Rock Cover - Unembedded	Flat Rock Cover - Total
Medium Stess Class
PT004WM	Small	+	0	0	++	0	+	0	0	0	0	++	++	--	+	+	0	0	+	0
PT002WM	Intermediate	+	+	--	+	0	-	-	0	+	+	++	++	--	+	+	0	-	+	0
PT001WM	Intermediate	0	0	0	0	+	-	--	0	0	+	++	++	--	+	+	0	-	0	-
+ + = Increasing Significant | + = Increasing Trend | 0 = No Change | - = Decreasing Trend | - - = Decreasing Significant

Results - Habitat - Biota Relationships

Correlations

Heatmap of spearman correlation coefficients of habitat and fish community endpoints, with green and blue indicating positive and negative relationships respectively. Strong Statistically significant correlations (rho > 0.6 & p < 0.1) are shown with red outlines. Moderate statistically significant correlations (rho > 0.4 & p < 0.1) are shown with orange outlines. Endpoints are arranged according to a Ward clustering.

Scatterplots

Scatterplots of habitat and fish community endpoints where absolute spearman correlation coefficients are greater than 0.3.

Site Specific Evaluations

By Les Stanfield November 10, 2025

As a follow up to the holistic analysis conducted by Patrick Schaefer, I asked for a subset of sites for which field crews had identified visually as becoming “more degraded over time. The intention was to provide a quality control check on the broader analysis. That is, are the broad patterns in the analysis in agreement with field impressions. Such intuitive interpretation provides a form of quality control for the conclusions.

A Crew leader (Mark Szonda) and Jeff Vandenburg (program supervisor) developed the following list of sites deemed to be undergoing “significant morphological change over the observed time of sampling.

The selected sites were:

• CC002WM
• DN008WM
• HU021WM (excluded due to lack of data)
• HU032WM
• DF005WM
• EC002WM
• HU024WM
• PT001WM
• DF008WM
• EC012WM
• HU029WM
• RG003WM

Interpretations from table changetbl1.xls were summarized to provide a synopsis of directions for each habitat endpoint. Each endpoint was classified as to whether over the 25-year period there was no change (0); nonsignificant increasing (+1) or decreasing (-1) change; or significant increasing (+2) or decreasing (-2) change. Habitat variables were ranked according to the mean change over time for all the sites for that category. For example, a value close to 2 indicates that nearly all the sites had significant changes in that variable. The results of this analysis is summarized below:

Ranked Habitat variables by amount of change over time

Habitat Variable	Amount of classified change
Non-Filamentous Algae Percent	1.82
Moss Percent	-1.64
Filamentous Algae Percent	1.36
D50	1.09
Wood Cover - Unembedded	1.00
Habitat Homogeneity	1.00
Flat Rock Cover - Unembedded	0.64
Wood Cover - Total	0.64
Usable Area	0.55
Depth - Mean	0.55
D50Max	-0.55
Bank Stability - Mean	0.45
Round Rock Cover - Total	-0.36
Flat Rock Cover - Total	0.27
Width:Depth Ratio	-0.27
Round Rock Cover - Unembedded	0.18
D85	0.18
Sediment Sorting Index	-0.18
Width - Mean	0.09
Bank Cover - Total	0.00
Bank Cover - Unembedded	0.00
Macrophyte Cover - Total	0.00

This analysis supports the broader conclusions that the most consistent measure of change over time in habitat is reflected in an increase in the amount of algae, both non and filamentous with a concurrent decrease in the amount of moss at a site. These conclusions are strong and consistent with ecological theory. Degrading streams have less stable substrates to which moss would attach. Algae is a much faster and seasonal plant and would be more adapted to grow on less stable substrates. Additionally, degrading streams are often found to receive increased nutrients which would also encourage algal growth.

While not significant, there were increasing trends in the mean point particle (D50), amount of unembedded and total wood cover and sites became more similar in depth and velocity profiles (Habitat homogeneity. This tendency is more complicated to interpret. A stream that is undergoing increased flow flashiness would be hypothesized to have greater movement of substrate and potentially loss of finer materials. This scenario would also result in increased erosion and therefore increased supply of wood into the channel. Finally, increased erosion typically manifests as either downcutting (degradation) or channel widening (aggradation). Assuming channel capacity doesn’t change, downcutting would tend to narrow the channel, decreasing mean width (as evidenced EC002WM). When stream beds are less erodible, the channel tends to widen as is evidenced at site CC02WM). There is clearly evidence of both scenarios in the dataset that complicate interpretation of broad patterns.

None of the geomorphic metrics consistently indicate a degradation in channel stability and in some cases the data are confounding. For example, in contrast to expectations, bank stability tended to increase at most of the sites over time. A cursory look at the 2025 site features data for these sites determined that at only one site (EC012WM) was there evidence of channel hardening, negating this activity as an explanation of these findings.

There were also no consistent changes in the width or depth, substrate or sediment sorting. At many of the site’s crews recorded evidence of sediment loading that could be a factor in their classification of these sites as being degrading over time. Perhaps there is a more nuanced and multivariate relationship between these metrics and channel stability, but there is no simplistic correlation between what biologists interpret as channel degradation and the OSAP derived metrics of channel stability. We will continue to seek an independent measure of channel degradation and to reevaluate the existing metrics of channel stability in attempts to improve our ability to measure this important stream metric.

The OSAP habitat module was designed to quantify conditions at a site over time and to provide guidance as to probable causal pathways. The strength of this module is hypothesized to be in its ability to help disentangle the effects of geomorphic processes at an individual site level, for example to understand the often opposite effects of degradation and aggradation. The two sites identified above will be used as case studies to ascertain whether observed patterns agree with ecological theory.

Site CC002WM: This site is on an intermediate sized stream significantly increased in both width and depth over the study period. This resulted in a predictable significant increase in usable area. There were no consistent changes in the traditional substrate metrics, which is supported by their being no change in rock cover, contrary to what I might have predicted. In contrast to expectations the bank stability indice does not decrease and in fact shows a positive tendency. A review of the cross-sectional profiles for this site shows a generally wide channel with shallow banks, many of which would be classified as point bars and therefore stable. The channel, while widening has maintained its overall form from year to year. I hypothesize that if calculated and compared there would also be an increase in the bankfull channel capacity over time. Like other sites there were significant increases in both filamentous and non-filamentous algae and a decrease in moss. There was also a significant increase in macrophytes. These increases imply considerable increase in nutrients and sunlight to the channel. I hypothesis that there would be a significant decline in the Hilsenhoff Index of benthic invertebrates for this site.

In general, the habitat data has documented that this channel is aggrading, but the channel has been able to accommodate this change with minimal impact on habitat quality. However, there is considerable evidence of eutrophication occurring at the site that I hypothesize is enhanced by a declining riparian cover.

Site EC002WM: This site is in the lower sections of Etobicoke Creek in an area that has cut through shale. The streams mean width and usable area have both significantly decreased over time suggesting serious degradation of the channel. The depth has only tended to increase but this could be a limitation of the protocol because in the later years points along several transects were “too deep to be safely measured” and therefore changes in depth were underestimated. Traditional substrate metrics change inconsistently. For example, the point particle size significantly increases over time as fewer smaller sands and gravels are observed. The stream is predictably becoming more poorly sorted, although not in a significant pattern. With downcutting the stream is predictably becoming significantly more homogeneous and the site is trending towards having more rock cover of both types.

Conclusions:

That the strongest signature across the dataset was in the significant decrease in moss that was coupled with significant increases in filamentous and non-filamentous algae was a surprise. This finding generates several other hypotheses that could be explored. For example, exploring whether other causal pathways explain this finding. As pointed out from the analysis above, there is value in interpreting the ecological changes in sites over time on both a broad and individual basis. There may be stream temperature or riparian vegetation changes that could explain the changes in vegetation.

It is also recommended that individual sites be evaluated for concerns on an annual basis.Such analysis could look for patterns, both on an annual basis and over time. The analysis should consider a holistic approach utilizing many of the main OSAP modules and other supporting data, such that ecological concerns could be identified and addressed in a timely manner.For example, analysis could evaluate benthic communities, temperature and habitat to see if flags emerge in any ecological factor (i.e., eutrophication, aggrading channels etc.). Where concerns are identified, methods could be modified to provide more diagnostic capabilities (greater numbers of points or more transects) or inclusion of other sampling methods in a truly adaptive manner. Conversely, sites which show no trends of change could be reconsidered for less sampling periodicity, or less intensity (i.e., rapid assessment modules). One recommendation then would be to generate an annual synopsis of site conditions that generates hypothesis that direct staff to action.

There is clearly more to the finding of the changes in vegetation that should be investigated.Leading hypothesis are that it truly is a measure of the stability of the bed, in that a rolling stone gathers no moss. However, equally plausible findings could be related to increased solar radiation (changes to riparian vegetation) and/or changes in stream temperature

Results and Conclusions

Model Fit

Overall models fit observed data well. Generally, the best fits were observed for Morphology endpoints, and weakest for some Substrate Sizes, Percent, and Cover endpoints. Most modeled endpoints contained some statistical deviations from expectations based on the selected distributional family, but visual inspection of residual plots showed good matches between observed and expected values to their theoretical quantiles.

Spatial Patterns

Spatial Correlations

At a high level, there are groups of endpoints with positive and negative correlations. Flat and round rock cover, substrate sizes, and algae and moss cover endpoints all showed some degree of positive correlations. Most of these endpoints were negatively correlated with wood cover and sediment sorting index. These endpoints were also positively correlated with reach slope, upstream catchment developed (%) and upstream catchment stress index, but uncorrelated with upstream agriculture. Stream morphology endpoints were also positively correlated with this cluster, but only weakly so.

Morphology

Overall there were clear and expected patterns in morphology across the Upstream Catchment Area (UCA) classes. For stream width, there was a non-significant pattern of low stress streams being narrower, especially among smaller streams. Conversely, for stream depth, there was a non-significant pattern of low stress streams being shallower, but only among small and intermediate streams. Together these patterns showed higher usable areas among small.

Substrate Sizes

Substrate sizes varied predictably with UCA, as smaller streams had smaller substrates. Substrates were also larger at high stress streams. Many of the monitoring locations with the smallest substrate sizes were found in Carruthers Creek watershed, but also some among the smaller streams in the upper Don and Humber watersheds.

Indices

Sediment sorting index was lowest in high stress streams, especially among the largest UCA. Mean bank stability was lowest in high stress streams, and highest in low stress streams. Within the low stress streams, there is a pattern of decreasing bank stability among larger streams than smaller streams. Habitat homogeneity varied by stream size, with larger streams having higher homogeneity. Upstream stress did not have a significant impact on homogeity.

Percent

Among the percent endpoints, there is a greater amount of variation among larger streams, than small and intermediate streams. Among small streams in particular, upstream stress seemed to play a large role in formation of algae which accumulated more in high stress sites then low stress sites. Filamentous and non-filamentous algae tended to be high in sites closer to lake Ontario and increasing further from the lake. Carruthers Creek and Frenchmans Bay watersheds locations were an exception and had high percent algae. Percent moss did not show this spatial pattern. High algae percentages may be related to elevated nutrient concentrations among more agricultural streams, especially in the upper Humber river watershed.

Cover

There is little difference in the spatial patterns between the unembedded and total fractions of wood, round rock and flat rock cover endpoints. There are also few consistent patterns among the cover endpoints in relation to overall UCA, but did show some weak patterns related to upstream stress. Wood cover tended to increase with decreasing stress, while flat and round rock cover increased with increasing stress. Wood cover tended to be lowest closer to lake Ontario, and higher further from the lake. This may be reflective of reduced forest cover in more densely urbanized areas closer to the lake. Round and flat rock covers both tend to be higher closer to lake Ontario, suggesting greater proportions of these cover types in more developed areas. These spatial patterns are confirmed in trend analysis which showed positive correlations with Upstream Catchment Stress Index and Upstream Catchment Developed (%).

Temporal Patterns

Temporal Correlations

Temporal trends in flat and round rock cover, substrate sizes and stream morphology were generally positively correlated, indicating share increasing or decreasing trends over time. Wood cover, bank stability, sediment sorting index and moss percent were mostly negatively correlated with the previous endpoints, indicating contrasting trends (i.e., one endpoint is increasing while the other is decreasing over time, or vice versa). Wood cover, bank stability, sediment sorting index and moss percent were generally positively correlated with each other.

Morphology

Stream width, depth and usable area are showing overall patterns of increasing over time throughout the TRCA. The rate of increasing width and depth was greatest in high and medium stress streams. Width to depth ratios were not changing consistently over time.

Substrate Sizes

The D50 endpoint showed consistently increasing trends over time, especially at large and intermediate streams. There are also patterns of increasing D50 at high and low stress classes, but not medium stress. Variation in temporal trends among D85 and D50Max were inconsistent across UCA and stress, but were both significantly increasing at high stress, intermediate sized streams.

Indices

Habitat homogeneity showed consistent increasing trends over time. The rate of increase was highest in high stress streams. Bank stability and sediment sorting index did not show consistent or significant temporal changes across UCA, stress or habitat stability classes. Spatially, bank stability tended to decrease among sites along, or direct tributaries to the Don and Huber rivers.

Percent

Filamentous and non-filamentous algae percents both increased over time at high and medium stress sites, and showed no consistent trends at low stress sites. The trends did not vary consistently by upstream catchment area. Rates of change were consistently higher for non-filamentous algae than for filamentous algae. Moss percent had an inverse pattern from algae, showing significant decreases over time at high and medium stress sites. Most areas of increasing algae were closer to lake Ontario, while the upper watersheds showed fewer patterns or even decreases.

Cover

There is little difference in the temporal patterns between the unembedded and total fractions of wood, round rock and flat rock cover endpoints. Round rock cover showed slight increasing trends, but only among high stress streams. Flat rock cover showed significant increases in small and intermediate, high stress streams. Spatially, sites showing increases in round and flat rock covers tended to be closer lake Ontario, while sites increasing in wood cover tended to be in the upper parts of the watershed.

Relationships Between Habitat and Fish Communities

There are no strong relationships (spearman rho > 0.6) between habitat and fish community endpoints. Some weak, but noteworthy relationships do exist however. Specifically, lithophilic spawners and open substratum spawners were positively correlated with a cluster of habitat endpoints including round and flat rock cover, substrate sizes, moss percent and habitat homogeneity. Sensitive darters were positively correlated with width, usable area, and width:depth ration, but not depth alone. Interestingly total community biomass, creek chub and community density were negatively correlated with these same morphology endpoints. Wood cover was also positively correlated with brook trout, creek chub, and coldwater stenotherm fish endpoints. Fish community richness, tolerant taxa richness and SATI biomass (a thermally weighted biomass measure) were positively correlated with width, depth and usable area.

Conclusions

Overall, morphology and percent moss and algae endpoints showed strongest and most consistent spatial and temporal patterns to broad landscape level drivers. Morphology also showed some of the strongest correlations with fish community endpoints. Cover endpoints showed less consistent spatial and temporal patterns, but did show some relationships with fish community endpoints.

The stress classes used in the current models are perhaps suboptimal for describing spatial patterns in stream habitat data. Though correlation between most habitat endpoints were stronger with the Upstream Catchment Stress Index than with Upstream Catchment Developed (%), perhaps different thresholds should be considered to better understand spatial patterns among in stream habitat endpoints.